{
 "cells": [
  {
   "cell_type": "code",
   "execution_count": 1,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "%matplotlib inline"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "\n",
    "# 1D Unbalanced optimal transport\n",
    "\n",
    "\n",
    "This example illustrates the computation of Unbalanced Optimal transport\n",
    "using a Kullback-Leibler relaxation.\n",
    "\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 2,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "# Author: Hicham Janati <hicham.janati@inria.fr>\n",
    "#\n",
    "# License: MIT License\n",
    "\n",
    "import numpy as np\n",
    "import matplotlib.pylab as pl\n",
    "import ot\n",
    "import ot.plot\n",
    "from ot.datasets import make_1D_gauss as gauss"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Generate data\n",
    "-------------\n",
    "\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 3,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "#%% parameters\n",
    "\n",
    "n = 100  # nb bins\n",
    "\n",
    "# bin positions\n",
    "x = np.arange(n, dtype=np.float64)\n",
    "\n",
    "# Gaussian distributions\n",
    "a = gauss(n, m=20, s=5)  # m= mean, s= std\n",
    "b = gauss(n, m=60, s=10)\n",
    "\n",
    "# make distributions unbalanced\n",
    "b *= 5.\n",
    "\n",
    "# loss matrix\n",
    "M = ot.dist(x.reshape((n, 1)), x.reshape((n, 1)))\n",
    "M /= M.max()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Plot distributions and loss matrix\n",
    "----------------------------------\n",
    "\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 4,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAZQAAADFCAYAAABzYARGAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADl0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uIDIuMi4yLCBodHRwOi8vbWF0cGxvdGxpYi5vcmcvhp/UCwAAIABJREFUeJzt3XmcjXX7wPHP11hm7EIqe5Fs06ghS3mqhygiPZSQsSWJpBQ9LUTKk6eiSDzKUgkRya9URGiRwWCMJVuMJdOUfZvl+/vjOjONMWbOmDPnPsv1fr3Oy5xz7nPua+45znXf3+X6GmstSimlVF4VcDoApZRSgUETilJKKY/QhKKUUsojNKEopZTyCE0oSimlPEITilJKKY/QhKKUUsojNKEopZTyCE0oSimlPKKg0wFkVq5cOVutWjWnw1BKKeWybt26P6y15XPazucSSrVq1YiOjnY6DKWUUi7GmN/c2U6bvJRSSnmEJhSllFIe4VZCMca0NsZsN8bsNMYMy+L5p4wxccaYTcaYZcaYqhmeizLG/Oq6RXkyeKWUUr4jxz4UY0wIMBFoCcQDa40xi6y1cRk22wBEWmtPG2MeA14HHjTGXAEMByIBC6xzvfav3ASZlJREfHw8Z8+ezc3LVIAJDQ2lUqVKFCpUyOlQlFJZcKdTvhGw01q7G8AYMxtoD6QnFGvt8gzb/wx0c/3cCvjWWvun67XfAq2BT3ITZHx8PCVKlKBatWoYY3LzUhUgrLUkJiYSHx9P9erVnQ7H9yUlwaZN8OOPsGMH1K8PTZpAnToQEuJ0dCpAuZNQKgL7M9yPB27JZvvewFfZvLZi5hcYY/oCfQGqVKly0RuePXtWk0mQM8ZQtmxZEhISnA7Ft504AcOGwfTpcPq0PBYWBmfOyM+lSsETT8ALL0Dhwo6FqQKTRzvljTHdkOatsbl5nbV2irU20lobWb581kOdNZko/QzkYPlyCA+HSZPgwQdh9mzYtw9OnYJff4UZM6BlSxg1Cho2hA0bnI5YBRh3EsoBoHKG+5Vcj13AGNMCeB5oZ609l5vXKqXywFoYMgTuvBMKFYJVq+CDDySpVK4MxkCNGtC9O3z6KSxaBAkJ0KgRvPGG09GrAOJOQlkL1DTGVDfGFAY6A4sybmCMaQBMRpLJkQxPfQ3cZYwpY4wpA9zleszvjB49mrp16xIeHk5ERARr1qxxOqSLFC9eHICDBw/SsWPHS2539OhR3n333Wzfq2nTpgCsWLGCtm3b5iqOhQsXEhf395iNl156iaVLl+bqPVQu/PvfkhgeewxiYqBZs+y3v/deiI2F9u0lEeXwWVDKbdbaHG/APcAOYBfwvOuxkUgCAVgK/A7EuG6LMry2F7DTdeuZ075uvvlmm1lcXNxFj3nTjz/+aBs3bmzPnj1rrbU2ISHBHjhwIM/vm5SUlOf3yKhYsWJubbdnzx5bt27dLJ/LHNPy5cttmzZtchVHVFSU/fTTT3P1Gnc5/VnwOW+8YS1Y26+ftampuXttUpK17dpZa4y1c+bkT3wqIADR1o1c4VbpFWvtl8CXmR57KcPPLbJ57QfAB+7sxx1PPiknYZ4UEQHjxl36+UOHDlGuXDmKFCkCQLly5dKfW7ZsGUOGDCE5OZmGDRsyadIkihQpkl5Cply5ckRHRzNkyBBWrFjBiBEj2LVrF7t376ZKlSp89NFHDB06lCVLllCgQAEeeeQRBg4cyLp163jqqac4efIk5cqVY/r06Vx99dUXxLVnzx66dOnCyZMnad++ffrje/fupW3btsTGxrJlyxZ69uzJ+fPnSU1NZf78+bz44ovs2rWLiIgIWrZsSZs2bXjxxRcpU6YM27ZtY8eOHRQvXpyTJ08CcPz4cdq0acPOnTu54447ePfddylQoMAF28ybN4/FixfTt29fFi1axPfff88rr7zC/PnzGTVqFG3btqVjx47ZHq+oqCi++OILkpKS+PTTT7nhhhs89ScOTDNnwtNPQ8eOMGGCNG3lRsGC0s9y113QrRuUKSN9LEpdJp0p74a77rqL/fv3c/3119O/f3++//57QEaf9ejRgzlz5rB582aSk5OZNGlSju8XFxfH0qVL+eSTT5gyZQp79+4lJiaGTZs20bVrV5KSkhg4cCDz5s1j3bp19OrVi+eff/6i9xk0aBCPPfYYmzdvvijZpHnvvfcYNGgQMTExREdHU6lSJcaMGcN1111HTEwMY8fK+In169czfvx4duzYcdF7/PLLL7zzzjvExcWxa9cuPvvss0v+bk2bNqVdu3aMHTuWmJgYrrvuuvTncjpe5cqVY/369Tz22GP897//zfE4BrVffoFevaTf5KOPLn8ocFgYfPEF3HAD3H8/7Nnj2ThVUPG54pA5ye5KIr8UL16cdevWsWrVKpYvX86DDz7ImDFjaNCgAdWrV+f6668HICoqiokTJ/Lkk09m+37t2rUjLCwMgKVLl9KvXz8KFpQ/xRVXXEFsbCyxsbG0dJ0tpqSkZJkwfvjhB+bPnw/Aww8/zNChQy/apkmTJowePZr4+Hjuv/9+atasmWVMjRo1uuT8jkaNGnHttdcC8NBDD7F69eps+2guZfv27dker/vvvx+Am2++OdukFfTOnoWoKLjmGvjsM3BdOV+20qVh8WKZq9KrFyxbBgX0XFPlnt8lFKeEhIRw++23c/vtt1O/fn1mzJhBgwYNLrl9wYIFSU1NBbhohn+xYsWy3Ze1lrp16/LTTz/lGFdOQ2m7dOnCLbfcwv/93/9xzz33MHny5PTk4G5MmfeRdj/j456oYpDWpBgSEkJycnKe3y9gvfQSbNsG33wj80o8oUoVeOst6N1bOukHDPDM+6qgoqchbti+fTu//vpr+v2YmBiqVq1KrVq12Lt3Lzt37gTgww8/5B//+AcgZfjXrVsHkH4VkZWWLVsyefLk9C/QP//8k1q1apGQkJCeUJKSktiyZctFr23WrBmzZ88G4OOPP87y/Xfv3s21117LE088Qfv27dm0aRMlSpTgxIkTbv/+v/zyC3v27CE1NZU5c+Zw6623AlChQgW2bt1KamoqCxYsSN/+Uu+f3fFSbvrxR/jvf+HRRz3f39GzJ9x9NwwdCq6/kVK5oQnFDSdPniQqKoo6deoQHh5OXFwcI0aMIDQ0lGnTptGpUyfq169PgQIF6NevHwDDhw9n0KBBREZGEpJN+3afPn2oUqUK4eHh3HjjjcyaNYvChQszb948hg4dyo033khERAQ//vjjRa8dP348EydOpH79+hw4kPX0nrlz51KvXj0iIiKIjY2le/fulC1blmbNmlGvXj2eeeaZHH//hg0bMmDAAGrXrk316tXp0KEDAGPGjKFt27Y0bdr0gia5zp07M3bsWBo0aMCuXbvSH8/ueCk3nD4NPXrI1cTYXM0ddo8x8L//yVyWHj0gJcXz+1ABzciIMN8RGRlpMy+wtXXrVmrXru1QRMqXBPVn4cUX4ZVXpI/jzjvzbz8zZ0ofzeTJ0Ldv/u1H+Q1jzDprbWRO2+kVilL+4OBBmbzYuXP+JhOAhx+Gpk1h+HAp26KUmzShKOUPhg+H5GQYPTr/92WMNKkdPiwd9Uq5SROKUr4uLk5qc/XvD1mM0MsXTZtChw7wn//AkSM5b68UmlCU8n3DhkHx4lJy3ptee03K3o8a5d39Kr+lCUUpX7ZypcxkHzYMMpT88YpateCRR+C993QYsXKLJhSlfNnzz8uM+EGDnNn/8OEyE3/ECGf2r/yKJhQ3JCYmEhERQUREBFdddRUVK1ZMv3/+/Pl82ef69etZsmSJW9veeuutxLgqZrZq1SrbSYtvvvlmtrPae/bsyfbt20lOTqZ06dJ5innBggXptcLUZfjhB1i9WiYaFi3qTAxXXQX9+kkRSa3zpXKgCcUNZcuWJSYmhpiYGPr168fgwYPT7xd2YxnVlMuYIJabhJLR119/TYkSJS75fHYJJSUlhWnTplGrVq1c7xcujrlDhw5uTZxUl/Cf/0DZslIOxUmDB0ttL12MS+XA/2p5OVG/Phv33nsvBw8e5OzZswwePJg+ffqQnJxMuXLl6NGjB9999x2TJ08mISGBZ555huLFi9O0aVP279/PwoULOXnyJAMGDCAuLo6kpCRGjhxJixYtGDlyJGfOnGHFihW88MILFxRjPH36NFFRUcTGxlKnTp0LEkSlSpWIjY0lJCSEBx54gIMHD5KSksKIESPYv38/R44c4bbbbqNChQosWbLkojifeeYZJkyYQL169QB44oknWLZsGddccw2zZ8+mbNmy3HrrrUyYMIGIiAgOHz7MrbfeyubNmy+K+ejRo8TGxjJu3Dj27NlDr169SExMpEKFCkybNo1KlSrRrVs3ypYty9q1azl8+DBvvPFG+kz8oBYbK30nL78MOdR+y3cVK8pqj++/L3XErrzS2XiUz9IrlDyaMWMG69atY+3atbz55pv89ddfABw7dozmzZuzadMmbrzxRvr3788333xDdHQ0hw8fTn/9yJEjad26Nb/88gvfffcdTz/9NMYYXnrpJbp27UpMTMxFlX0nTJhAmTJl2Lp1Ky+88AIbslgb/Msvv6RatWps3LgxvXLx4MGDufLKK1m1alX6CooZ42zSpMkF73Hs2DGaNWvGli1baNKkCaOyGe0TFhaWbcz9+/enT58+bNq0iU6dOl1QkfnIkSP88MMPLFy4kOeee87NIx/gXn9dEsnjjzsdiXjmGTh3Dt5+2+lIlA/zvysUJ+rXZ+Ott95i0SJZETk+Pj594arChQunn2nHxcVRq1YtqlatCkgJ+JkzZwLwzTff8NVXXzFmzBhAqvbu27cv232uXLmSZ599FoAGDRpQt27di7YJDw9n2LBhDBs2jHvvvZdml1gWNmOcmRUsWJBOnToB0K1bN7p06ZJtXNlZs2YNixcvBqB79+68+OKL6c/dd999GGMIDw+/ZE2yoPLbbzBrFjzxhDR5+YJatWS9lIkT4dlnoWRJpyNSPkivUPJg6dKlrFy5kp9//pmNGzcSHh6e3vwUFhaWY2l5kFL1CxcuTO+T2bdvX/p6IXlRu3ZtoqOjqVu3LsOGDePVV1/Ncjt344S/y9VnV5r/chTJsJ6Hr9WWc8Qbb0ifxVNPOR3JhYYOhaNHYcoUpyNRPkoTSh4cO3aMK664grCwMLZs2cLatWuz3K5OnTps376d/fv3Y61lzpw56c+1atWKd955J/1+WvNVdiXmmzdvzqxZswDYuHFjlqXtDxw4QPHixXn44Yd5+umnWb9+fY7vm1lycnL6QlezZs1KL1ufsTT/vHnz0rfP7r0bN27M3LlzAfjoo49o3ry5WzEEnb/+gqlToWtXqFTJ6Wgu1LCh1BF76y1ISnI6GuWDNKHkQZs2bTh9+jR16tThhRde4JZbbslyu6JFizJhwgRatGhBZGQkpUuXppRrYaThw4dz6tQp6tevT926dRnhGu9/5513snHjRho0aHDBlzbAgAEDSExMpHbt2owaNSrLhb42btxIw4YNiYiI4NVXX+Xf//43AH379qVFixa0aNEix9+vVKlSrFq1irp167J69WpecM3UfuaZZxg/fjw33XRTep9RTjFPnDiRKVOmEB4ezpw5c3hLa0Rlbdo0mZ2ew6qfjnnySSlU+fnnTkeifJCWr/eSkydPUrx4cay1PProo9SvX5+BAwc6HZbfCYTPwiWlpkLNmjKRcdUqp6PJWkoK1KgBVavCihVOR6O8RMvX+5hJkyYRERFBnTp1OHPmDI888ojTISlfs2QJ7N7t28vvhoRIkcrvv4fNm52ORvkYvUJRfiWgPwtt2sCGDbB3L7gxYdYxiYnSv9OjB0ya5HQ0ygsC7grF1xKf8r6A/gzs3AlffSUrJPpyMgEZyvzQQ7Ky49GjTkejfIhfJJTQ0FASExMD+wtFZctaS2JiIqGhoU6Hkj8mTZLmJH9ZcnfAAFnjfsYMpyNRPsQvmrySkpKIj4/3yJwH5b9CQ0OpVKkShQoVcjoUzzp9WsqbtGolRRj9RdOm8McfsG2bzJtRAcvdJi+/mClfqFAhqlev7nQYSuWPOXOk6ah/f6cjyZ3HH4du3eC778CNYegq8OlphVJOmzoVbrgBbrvN6Uhy51//gjJlpGikUriZUIwxrY0x240xO40xw7J4vrkxZr0xJtkY0zHTcynGmBjXbZGnAlcqIMTFwY8/Qp8+4GYJHJ8RGgoPPwyffSYjv1TQyzGhGGNCgInA3UAd4CFjTJ1Mm+0DegCzsniLM9baCNetXR7jVSqwvP8+FCokX8z+qHdvOH8ePvrI6UiUD3DnCqURsNNau9taex6YDbTPuIG1dq+1dhOQmg8xKhWYzp2TUVLt2/vvGiPh4dCoEfzvf+BjA3yU97mTUCoC+zPcj3c95q5QY0y0MeZnY8x9WW1gjOnr2iY6ISEhF2+tlB/7/HNpKurTx+lI8qZPH9iyBdascToS5TBvdMpXdQ036wKMM8Zcl3kDa+0Ua22ktTayfPnyXghJKR8wdSpUqeL/I6Q6d5bFwKZOdToS5TB3EsoBoHKG+5Vcj7nFWnvA9e9uYAVwcWlcpYLNnj3w7bfQq5dMaPRnJUrAgw/KHBo3l0ZQgcmdhLIWqGmMqW6MKQx0BtwarWWMKWOMKeL6uRzQDIi73GCVChjTpsmorp49nY7EM/r0gVOnZE6NClo5JhRrbTIwAPga2ArMtdZuMcaMNMa0AzDGNDTGxAOdgMnGmLQVn2oD0caYjcByYIy1VhOKCm6pqdIZ37KlNHkFgsaNZS7N9OlOR6Ic5NZMeWvtl8CXmR57KcPPa5GmsMyv+xGon8cYlQos338P+/bBmDFOR+I5xkj14WHDpNBljRpOR6QcoDPllfK26dOhZEm4L8tBj/6rWzep6aUFI4OWJhSlvOnECZg3Tzqxw8KcjsazKlaUZryZM6VZTwUdTShKedP8+VJduEcPpyPJHz16SHOeLg8clDShKOVN06fLuvFNmjgdSf5o3x5KldLO+SClCUUpb9mzRzrke/Twv0KQ7goLk4mO8+frnJQgpAlFKW+ZOVMSib8WgnRXVJQ0682b53Qkyss0oSjlDdbK6Kc774TKlXPe3p81bgzXX6/NXkFIE4pS3vDDD9LkFRXldCT5zxjo3h1WroS9e52ORnmRJhSlvGHmTCmg2KGD05F4R7du8q+ukxJUNKEold/OnoW5c+H++6F4caej8Y6qVeEf/5BEquukBA1NKErlty++gGPHpBkomHTvDr/+Cr/84nQkyks0oSiV32bOlFnkd9zhdCTe1bGjrDs/c6bTkSgv0YSiVH46cgS++gq6dvX/dU9yK61e2ezZstyxCniaUJTKT598AikpgT/35FK6d4c//4Qvv8x5W+X3NKEolZ8+/BBuugnq1XM6Eme0bAkVKshxUAFPE4pS+SUuDtatC96rE4CCBaFLF1i8GBITnY5G5TNNKErll5kzpd/koYecjsRZ3btDUpIuDxwENKEolR9SUuDjj6F1a2nyCWY33gj162uzVxDQhKJUflixAuLjg2/uSVbSSrH8/DPs2OF0NCofaUJRKj/MnCnDZu+91+lIfEOXLrI8sF6lBDRNKEp52qlTsh7IAw8E3jK/l+uaa6BFC6ntpcsDByxNKEp52oIFklS0uetC3btL9eHVq52OROUTTShKedrMmVCtGjRr5nQkvuW++6TispZiCViaUJTypAMHYNkymXtSQP97XaBYManv9emncOaM09GofKCfeKU86eOPpY8gmCczZqd7dzh+HD7/3OlIVD7QhKKUp6Qt89u0KdSs6XQ0vun222UJ5BkznI5E5QNNKEp5yrp1Um4lGJb5vVwFCshVyjffwMGDTkejPEwTilKeMn26rP/xwANOR+LboqKkWVCXBw44biUUY0xrY8x2Y8xOY8ywLJ5vboxZb4xJNsZ0zPRclDHmV9dNT91UYDp3TkrV33cflC7tdDS+rWZNaRacMUOXBw4wOSYUY0wIMBG4G6gDPGSMqZNps31AD2BWptdeAQwHbgEaAcONMWXyHrZSPmbxYln3Q5u73BMVJc2D0dFOR6I8yJ0rlEbATmvtbmvteWA20D7jBtbavdbaTUDmKbCtgG+ttX9aa/8CvgVaeyBupXzLjBkyG7xlS6cj8Q8PPCDNg9o5H1DcSSgVgf0Z7se7HnOHW681xvQ1xkQbY6ITEhLcfGulfMTvv8uKhN26Bd8yv5erdGlpHpw1S5cHDiA+0SlvrZ1irY201kaWL1/e6XCUyp1Zs6RcvTZ35U6PHvDXX/DFF05HojzEnYRyAKic4X4l12PuyMtrlfJ91sIHH0DDhlAnc9eiylaLFtJMOG2a05EoD3EnoawFahpjqhtjCgOdgUVuvv/XwF3GmDKuzvi7XI8pFRjWroXYWOjd2+lI/E9IiFylLFkiJWuU38sxoVhrk4EBSCLYCsy11m4xxow0xrQDMMY0NMbEA52AycaYLa7X/gmMQpLSWmCk6zGlAsP770uJ+s6dnY7EP/XqJXNSpk93OhLlAcb62DjwyMhIG61DCZU/OHUKrr4aOnTQ0Up5cccdsG8f/PqrFtT0UcaYddbayJy207+eUpdr3jw4cUKbu/Kqd2/YvRu+/97pSFQeaUJR6nK9/77M+r7tNqcj8W//+heUKiXHU/k1TShKXY4dO2DVKukDMMbpaPxbWJisOT9/Phw96nQ0Kg80oSh1OT74QEYp6dwTz+jdG86elTk9ym9pQlEqt5KSpBP+nnukU17l3U03QUQETJ2qBSP9mCYUpXLr88/h8GF49FGnIwkcxkDfvrBhg8ztUX5JE4pSufXee1ClCrTWOqce1bWrrDv/3ntOR6IukyYUpXJjxw5YtkzOprUQpGeVLClJZfZsqfGl/I4mlABz+DCMGwePPy6rrCYnOx1RgJk8GQoW1Lkn+aVfPzhzBj780OlI1GXQhBIgli6FVq2gYkUYPFgGIbVqBZUrw9NPw7FjTkcYAM6ckRIhHTrAVVc5HU1gatAAbrlFmr20c97vaEIJAJ9/DnffDdu2wXPPyUJ4f/0lE7mbNIG335bnT5xwOlI/9+mnsipjv35ORxLY+vWDrVth5UqnI1G5pLW8/NxXX0H79jLq8ttvoUSJi7dZsAA6dYJmzWT7okW9H2dAaNoUEhMlc+tkxvxz5oyUtW/dGj75xOloFFrLKygsWwb33w/160sF8KySCUgLzccfw+rVknzOnvVunAEhJgZ++knOnjWZ5K+wMClrP3++dAoqv6EJxU/t2SPJoUYN6XwvXTr77R98UPpVli6F/v29E2NAGT9eLu169HA6kuDQv7+MKNEhxH5FE4ofshYee0xOlP/v/6BsWfdeFxUFw4bJAnnLl+dvjAHlyBEpCdKjB5Qp43Q0waFmTWjTBiZN0jXn/YgmFD80ezZ8/TWMHi3z63LjpZfguutkkrc2fbnpvffg/Hl44gmnIwkugwZJMp892+lIlJs0ofiZP/+EJ5+UJcwffzz3rw8Lk+/HX3+FV1/1fHwB59w5ePddGSZXq5bT0QSXf/4T6taViVU+NnhIZU0Tip8ZOlQGGk2ZcvkTtVu0gG7dYMwYGWKssjF3Lvz+u2Rx5V3GyHGPidEhxH5CE4of+eknKcY6eLAUZs2LN9+UUWGPPaYnf5dkLbz1FtSuDS1bOh1NcOraVToJx41zOhLlBk0ofuSFF+DKK2HEiLy/V/nyMGqUnPh9+23e3y8grV4t1W8HDdKhwk4JC5MOv88/l2WClU/ThOInVqyA776TmfDFinnmPXv3lk79l17Sq5QsjR0LV1wBDz/sdCTBrX9/qZ/25ptOR6JyoAnFD1gLw4fLWk6eXIKjSBF4/nlYs0Zm0KsMNm+GL76QqxMtLeCsihWhe3dZc/73352ORmVDE4ofWL5cmqaee05aADypRw+oVk0Sll6lZDBmDBQvDgMGOB2JAnj2WRlxp30pPk0Tio+zVpqkKlaERx7x/PsXLgwvvgjR0bB4seff3y/t2iVzH/r1kyYv5bzrr5eCdO++C0ePOh2NugRNKD7u22/hhx+kaSo0NH/28fDDMtlR+1Jcxo6VNvunnnI6EpXRc8/B8eOSVJRP0oTi4155RdY06dUr//ZRqJBcpcTEwJdf5t9+/MLBg1KbpmdP6bRSviMiQiaYjhsHp087HY3KgiYUH/bTT7BqlSyQVaRI/u6rSxdJXK+/nr/78XlvvSVFCZ991ulIVFb+/W9ISJAOeuVzNKH4sNdfl1qE3lhttlAhaeFZuRJ+/jn/9+eTDh+GiRPhoYfg2mudjkZl5dZboXlzeO01WTdF+RS3EooxprUxZrsxZqcxZlgWzxcxxsxxPb/GGFPN9Xg1Y8wZY0yM66a1qN20bZvM5Xr8cRls5A19+kgCC9qrlNGjpQikJ2aOqvzzyitw6JAkf+VTckwoxpgQYCJwN1AHeMgYUyfTZr2Bv6y1NYC3gP9keG6XtTbCddO1U930xhvSzDVwoPf2Wby4JLCFC2H7du/t1yfs3QuTJ8vlYI0aTkejsnPbbbKa42uvSSe98hnuXKE0AnZaa3dba88Ds4H2mbZpD8xw/TwP+KcxWqvich06BDNnSr/wlVd6d98DB8pQ4v/+17v7ddzIkVCggIxOUL7vlVek9LbOnvcp7iSUisD+DPfjXY9luY21Nhk4BqQt+1TdGLPBGPO9Mea2rHZgjOlrjIk2xkQnJCTk6hcIROPHS7/w0097f99XXimJbOZMSWxBYds2mDFDLs8qVXI6GuWOm2+Gjh3lUv6PP5yORrnkd6f8IaCKtbYB8BQwyxhTMvNG1top1tpIa21k+fLl8zkk33b8uCxS969/ydwQJwwZIglt/Hhn9u91L70k5VWGXdQ9qHzZyJEyfHjMGKcjUS7uJJQDQOUM9yu5HstyG2NMQaAUkGitPWetTQSw1q4DdgHX5zXoQDZliiQVJ0etXnedJLT33guCJuqff4ZPP5U1AYL8ZMbv1K4ts3InTNBKxD7CnYSyFqhpjKlujCkMdAYWZdpmERDl+rkj8J211hpjyrs69THGXAvUBPQvfwnnz8s0iDvvhMhIZ2N55hk4dgz+9z9n48hXqanSaXTNNTqB2FQAAAAQf0lEQVTvxF+NHi1VDYYMcToShRsJxdUnMgD4GtgKzLXWbjHGjDTGtHNt9j5Q1hizE2naSms7aA5sMsbEIJ31/ay1f3r6lwgUH38sE7V94butYUO44w5JcOfPOx1NPpk+XYqYvf6698ZmK8+qWFHqEi1YAEuXOh1N0DPWx4o3RUZG2ujoaKfD8LrUVKhXT0ZYbdjgG+s5ff21jM6cNk2qEgeUY8ek4GCNGrKQli8ccHV5zp6VtedDQ6V+UKFCTkcUcIwx66y1Obab6Ex5H/F//wdbt8rVia98t911F4SHS63E1FSno/GwkSOlhMfbb/vOAVeXJzRULqXj4rRwpMM0ofiI11+HqlWlQrevMEYSXFxcgBWNjIuTRNKnjww/Vf7v3nuhVStZ2EcX4XKMJhQfsHq13J56yveu1h94QBLdq68GSGn7lBSZDV+ypHToqsBgjIxzP3NGF0VzkCYUH/Dyy1Chgpww+5pChWDoUKl8vGyZ09F4wLhxMlT4nXd0mHCgqVVL/jPNmydDwZXXaae8w378EZo1k1InTsyMd8e5c9J3Xa2aVCP22y6HHTvgxhulc2jhQj/+RdQlJSdDkybw22+wZYueNHiIdsr7iZdfls98Px8um1mkiEwiX71a1rf3SykpskpZaKjM2NRkEpgKFpRhiUePereyqgI0oTjqp5/gm29kEmGxYk5Hk73evWX+34gRftqX8vbbspby+PG6EmOgq1dPinzOmSPNX8prNKE46OWXoVw56N/f6UhyFhoqVymrVsGKFU5Hk0u//CIdQe3bS6kOFfiGDZNyE336aFkWL9KE4pA1a2Ti4JAhvn91kuaRR+Tk3q+uUv78U4aqVawoTSHa1BUcChWCuXPl792pk0x+VPlOE4oDrJVEUr68VEz3F6GhsqT3ypUyEdPnpaZCVJTUs5k7V5ajVMGjenVZlmD9ehmTr/KdJhQHzJ8vHdyvvOJ/JaQefVRGZz79NCQlOR1NDsaOhcWLZRGmhg2djkY5oV076aScNAlmzXI6moCnCcXLzp6Vz3d4uHR0+5tChWRNox07fLzKxYIF8Nxz0tzlT5eByvNGj4Zbb5X/cD/95HQ0AU0TipeNGyfLl7/5JoSEOB3N5bnnHmjZUvpSEhOdjiYLP/0EXbpAo0bab6LkLOizz2Q1znvvlbMhlS80oXjR4cNSwqRdO/jnP52O5vIZIwnx+HEZqeZTduyQL41KleCLL2QlRqXKl4evvpIP7913w5EjTkcUkDSheNFzz0mpobFjnY4k7+rVg759pdlr82ano3E5cEC+LIyRLw+dJa0yqlFD+tQOHYI2bWQJA+VRmlC8ZPFiWc9pyBBZhiMQjBoFZctC9+4+sAjXb79B8+Zy5rl4sXx5KJXZLbfIhMeNG6WZwCfbbP2XJhQvSEiQ+VXh4dLvECjKlYMpU2RNI0ebvnbuhNtukzknS5fKl4ZSl3LvvdKnsnmzLEuqzV8eowkln1krdbr++gs+/FDqYgWS9u2hZ08YM0YKXXpdbKxcmZw5I4XGNJkod7RtK1eyO3fCP/4B+/Y5HVFA0ISSzz76SE6GRo2SK5RANG4cVKkiTV8nT3pxxwsWQOPG8vOKFRAR4cWdK7/XsiUsWSITXxs2lLpCKk80oeSjzZtlCsRtt/luaXpPKFlSJiTv3i3lWfJ9ueDUVGk7vP9+WUs8Olr+VSq3mjeXOkilS8Odd8LkyU5H5Nc0oeST/ftlwFGJEvDxx/4758RdzZvLkOjZs6U8S745dEjGXb/8spRV+f57KYOs1OW64QZJKi1bSvt0jx46AuwyaULJB0ePSjI5flxGr1au7HRE3jF0KDz2GPznPzBxooff3FopnVG3riwd+c47MmkxNNTDO1JBqXRpmbf04ovSTl2vnlRvVbmiCcXDzp6Vlpjt26WJP1D7TbJijHzPt2snaxstWOChN961Czp0gK5d5WwyJkbWDdcZ8MqTQkJg5EiptFCiBLRuLYuyHT7sdGR+QxOKB/3+uzTDLl8OH3zg37PhL1dICHzyiVQ96dRJavJdtsREGDwYateGb7+F11+XjtNatTwWr1IXadhQKhQPHSpDM2vUkCbWU6ecjsznaULxkM2b5Us0JkYWiQvmdZyKFpXv/7vvlsXDBg6Upb7ddvAgPP88XHedrLQYFQW//ipVNQO9M0r5htBQGQu/dat8kEeMkMQyZozMd1JZ0oSSR9ZKp3vTpvKluWoV/OtfTkflvBIlYOFCGd02YYIUlMx2qL+10jEaFQXVqsFrr8kl3qZN8L//ace7ckaNGvDppzLJql49qZ9UubKcJcXGOh2d77HW+tTt5ptvtv5i40Zrmze3Fqy95RZr4+Odjsg3TZ1qbViYtUWLWvvKK9aeOeN6IjXV2q1brR0+3NoaNeRAFitm7cCB1u7c6WTISmVt40Zre/SwtlAh+byGh1v7+uvW7t3rdGT5Coi2bnx/G+tja7lGRkba6Ohop8PIVkyMFEX84AMZHPLaa9J3p60xl/bbb3K1smL+H3S86gcGXv81tX/7igK/7ZXO9TvugG7dZERDqVJOh6tU9hISZBXQjz6Cn3+Wx2rXlo78Vq1kwm0AfY6NMeustZE5budOQjHGtAbGAyHAVGvtmEzPFwFmAjcDicCD1tq9rueeA3oDKcAT1tpsx+L5akLZt08m1U6dCmvXShNr794yKOSKK5yOzgdZK6NjYmOlgykmRkbP7NwJwEmKsSLknxxtfDfVB7bl5vaVdASw8k87d8qQ4yVLZF7UuXNyklSnjiSWG2+E+vWlyaxcOaejvSweSyjGmBBgB9ASiAfWAg9Za+MybNMfCLfW9jPGdAY6WGsfNMbUAT4BGgHXAEuB6621KZfan9MJJe17cOtW2LYNNmyQUVu7dsnzderIMrjdugVxIjl3TkZg/fGH3A4dktvBg7J62O7dsGePTMRJc9VV8p+rSRNs4yasLXAL731QmNmzpQxXkSLQrJksrFe7ttyuvx7Cwhz7LZXKvVOn5MQp7bZmzYWd+GXKwLXXyq1qVekbvPpquZUtKwmnbFlZFMyHeDKhNAFGWGtbue4/B2CtfS3DNl+7tvnJGFMQOAyUB4Zl3DbjdpfaX14Sytlj51jz9FykcRNSrVTpsKmQkiK3pGRIOi/roZ87J19mZ87CieNw7DgcOwrJGdJd0TCZ+lCnjtwqVfLg9IdLHfvMj6f9QpmfS3s88y019e9/024pKX//m5IiIwhSUuRAJCVJ/fnz5+WgnD0rt9On/76dOCEJ4vhxeS4rRYtKh3raf5gaNeTMrG7dS65NcuIErFwJ330n8xU3bbrwVyxTRnJRhQrSvFiypHT4FysmSSjtVrCgNDkWLAgFCvx9M+bvv1fGn9PuZ/WzUh5jLWHHDlM6PpbS8bGU+H0nxRN2U/zIboon/kZI0rksX5ZcuCjnw0qSFFaS5NDiJBcuSnKRYqQUDiOlUCgphUJJLViE1IKFSQ0pJP8WKIgNKYgtEIItEEKq619rQqj/Tl9Cr7j8xebcTSgF3XivisD+DPfjgcwlXdO3sdYmG2OOAWVdj/+c6bUVswi2L9AXoEqVKm6ElLUziaf5x/vdL/v1Wb8psMF1CzQFCsg3cOHCckZUuLC05aV9SxcrJkmiTBn5Fi9ZUm6lSl14NnX11XKmVaJErr+ZS5SQtY7atJH7Z87ICOFt2+TfQ4dkfs/hw9KykJbXTp+W3KeUbzPA1a5by0zPWUpzlGs4yNUcoiyJ6beS549T6vwxSh07RjFOUZTTFOUoxThAEc5RhHOEcpZCJFGY8xTmPCGkEELWhfQSnu+ap4TiLncSSr6z1k4BpoBcoVzu+5SqXJKDK3dizN9npyEhf9/Svjt9qvP8Ul/AmR/PfKqd+fG0W8bT8pCQvx8LCfn737TTeR88LQ8Lk+oC7lQYsPbvK820i67k5L8vzlJSLrywy3xxl9XPSnmPAcq4brkrbpoMZFnY2/XBN6kp6a0RxqZSuUbJPEfrDncSygEgYzWqSq7Hstom3tXkVQrpnHfntR5ToFAI19x2XX69vfIxxsgJQuHCTkeilK8o4Lo50wfjzsTGtUBNY0x1Y0xhoDOwKNM2i4Ao188dge9cY5cXAZ2NMUWMMdWBmsAvngldKaWUL8nxCsXVJzIA+BoZNvyBtXaLMWYkMtllEfA+8KExZifwJ5J0cG03F4hDrtIez26El1JKKf+lExuVUkply91RXlrLSymllEdoQlFKKeURPtfkZYxJAH7L49uUA/7wQDj+TI+B0OOgxyCNHgdxOcehqrU269nJGfhcQvEEY0y0O+19gUyPgdDjoMcgjR4HkZ/HQZu8lFJKeYQmFKWUUh4RqAllitMB+AA9BkKPgx6DNHocRL4dh4DsQ1FKKeV9gXqFopRSyss0oSillPKIgEooxpjWxpjtxpidxphhTsfjLcaYysaY5caYOGPMFmPMINfjVxhjvjXG/Or6t4zTseY3Y0yIMWaDMWax6351Y8wa12dijqvAaUAzxpQ2xswzxmwzxmw1xjQJts+CMWaw6/9CrDHmE2NMaDB8FowxHxhjjhhjYjM8luXf3oi3XcdjkzHmprzuP2ASimup4onA3UAd4CHXEsTBIBl42lpbB2gMPO763YcBy6y1NYFlrvuBbhCwNcP9/wBvWWtrAH8BvR2JyrvGA0ustTcANyLHI2g+C8aYisATQKS1th5S1LYzwfFZmA60zvTYpf72dyMV4GsiCxxOyuvOAyahIOvW77TW7rbWngdmA+0djskrrLWHrLXrXT+fQL5AKiK//wzXZjOA+5yJ0DuMMZWANsBU130D3AnMc20SDMegFNAcqQCOtfa8tfYoQfZZQCqph7nWZyoKHCIIPgvW2pVIxfeMLvW3bw/MtOJnoLQx5uq87D+QEkpWSxVftNxwoDPGVAMaAGuACtbaQ66nDgMVHArLW8YBz0L6OqhlgaPW2mTX/WD4TFQHEoBprqa/qcaYYgTRZ8FaewD4L7APSSTHgHUE32chzaX+9h7/zgykhBL0jDHFgfnAk9ba4xmfcy14FrBjxI0xbYEj1tp1TsfisILATcAka20D4BSZmreC4LNQBjn7rg5cAxTj4magoJTff/tASiheXW7Y1xhjCiHJ5GNr7Weuh39Pu4R1/XvEqfi8oBnQzhizF2nuvBPpSyjtavaA4PhMxAPx1to1rvvzkAQTTJ+FFsAea22CtTYJ+Az5fATbZyHNpf72Hv/ODKSE4s5SxQHJ1VfwPrDVWvtmhqcyLs0cBXzu7di8xVr7nLW2krW2GvK3/85a2xVYjixLDQF+DACstYeB/caYWq6H/omsmBo0nwWkqauxMaao6/9G2jEIqs9CBpf62y8CurtGezUGjmVoGrssATVT3hhzD9KOnrZU8WiHQ/IKY8ytwCpgM3/3H/wb6UeZC1RBlgR4wFqbucMu4BhjbgeGWGvbGmOuRa5YrgA2AN2steecjC+/GWMikIEJhYHdQE/k5DFoPgvGmJeBB5ERkBuAPkj/QEB/FowxnwC3IyXqfweGAwvJ4m/vSrYTkObA00BPa22elssNqISilFLKOYHU5KWUUspBmlCUUkp5hCYUpZRSHqEJRSmllEdoQlFKKeURmlCUUkp5hCYUpZRSHvH/mjQ12FO9YI8AAAAASUVORK5CYII=\n",
      "text/plain": [
       "<Figure size 460.8x216 with 1 Axes>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    },
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAWAAAAFgCAYAAACFYaNMAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADl0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uIDIuMi4yLCBodHRwOi8vbWF0cGxvdGxpYi5vcmcvhp/UCwAAIABJREFUeJzt3XmcHHWd//HXp3tyEBIIEAyBAANyLSSAIcrNAkFBwI0/UeRQoiL3fSyni6DLobIQlgU0giy3CAgIcqywIAgkmhCW+wwEEo4EDHckycz398enerpn0pPpmenqb1f3+/l49KNTXTVVHzrMp9+p/lZ9LYSAiIjUXi52ASIizUoNWEQkEjVgEZFI1IBFRCJRAxYRiUQNWEQkEjVgkSZhZtub2Qux65AiNWBpema2n5lNN7OPzewtM7vbzLbr5z5fM7NdqlVjBccLZrbesrYJITwcQtiwj/t/zcwWmdmILq/PTI7d2pf9Njs1YGlqZnY8MBk4BxgJrAVcCkyMWVe1mVlLFXbzKrBvyT7HAkOqsN+mpQYsTcvMVgR+AhwRQvh9COGTEMLiEMIdIYR/TbYZZGaTzezN5DHZzAYl60aY2Z1m9r6Z/d3MHjaznJldgzfyO5JUfVKZY+9oZnPM7CQzm5ck76+b2e5m9mKyv9NKtv+SmT2WHOstM/svMxuYrHso2ez/kuN9u2T/J5vZ28CVhdeSn/l8coxxyfLqZjbfzHZcxlt2DXBAyfIk4Oo+vfkCqAFLc9saGAzcuoxtTge2AjYHNgO+BPwoWXcCMAdYFU/PpwEhhPBd4HXgayGEoSGEn3ez79WS468BnAH8GvgOsAWwPfBvZrZOsm0bcBwwIql7AnA4fsAdkm02S453Y8n+VwbWBg4uPXAI4RXgZOBaMxsCXAlcFUJ4cBnvxVRgBTP7JzPLA/sA1y5je+mBGrA0s1WAd0MIS5axzf7AT0II80II84GzgO8m6xYDo4C1k+T8cOjdzVUWA2eHEBYDv8Wb60UhhI9CCM8Az+JNnxDCjBDC1BDCkhDCa8CvgH/uYf/twI9DCJ+FEBZ2XRlC+DXwMjAt+e84vYKaCyn4y8BzwNwKfka6oQYszew9YEQP50dXB2aXLM9OXgP4Bd7A/sfMZpnZKb09fgihLflzoUG+U7J+ITAUwMw2SE53vG1mH+LnrDt9IVbG/BDCP3rY5tfAGODiEMJnFdR8DbAf8D10+qHf1IClmT0GfAZ8fRnbvIn/E75greQ1kqR6QghhXeBfgOPNbEKyXbVvM3gZ8DywfghhBfx0h/XwM8uswcyG4l9AXgGcaWYr91RECGE2/mXc7sDvK6hblkENWJpWCOED/NzrJckXYEPMbICZfdXMCudtbwB+ZGarJkOwziA572lme5rZemZmwAf4edr25OfeAdatYrnDgA+Bj81sI+CwLuv7cryLgOkhhB8CfwR+WeHPHQjsHEL4pJfHky7UgKWphRD+Azge/2JtPvAGcCRwW7LJvwPTgSeBp4DHk9cA1gfuAz7G0/SlIYQHknXn4o37fTM7sQqlnoj/0/8j/LTBjV3WnwlclRxv7552ZmYTgd0oNvLjgXFmtn9PPxtCeCWEML0XtUs3TDdkFxGJQwlYRCQSNWARkUjUgEVEIlEDFhGJpBo36JAGMGLEiNDa2hq7DJFMmTFjxrshhFX7+vNqwAJAa2sr06drZJFIb5jZ7J636p5OQYiIRKIGLNJMPvgAHn0UZs8GXQMQnRqwSKMLAW66CcaMgeHDYdttobUVPvc5OOUU+Oij2BU2LTVgkUb29tswYQLsvTfk83D22XD77XDppbDTTvCzn8GGG8K998autCnpSziRRjVnjjffOXPgkkvgkEO8CRccdhhMmwYHHQRf+xrccAPstVe8epuQErBII5o7F7bf3hPw//wPHH545+ZbsOWW8PDD8MUvekq+6aba19rE1IBFGs3ixd5M330X7r/fz/kuy4or+imIrbeGSZPg6adrU6eoAYs0nJNO8pEOV1wB48dX9jNDh8LNN3sz3msv+PDDdGsUQA1YpLHceSdMngxHH+0puDdWWw1uvBFeeQWOOCKd+qQTNWCRRvHxx36ud8wY+MUv+raPHXaA006Da6+FP/2puvXJUtSARRrFmWfCG2/Ar34FAwf2fT+nnQbrr+/N/B89zekp/aEGLNIInnzSTz0cdBBss03/9jV4sI8TfvllOPfc6tQnZakBizSCk06CFVaA886rzv522QW+/W0/lfHmm9XZpyxFDVgk6x54wIeRnXYarNzjzPKVO+ccWLIEfvKT6u1TOlEDFsmyEPx+DqNHw5FHVnff667rV89dfjm8+GJ19y2AGrBItt16K/z1r3DWWX7uttp+9CPf77/9W/X3LWrAIpkVAvz7v/uIhQMOSOcYI0fCscf6JcrPPZfOMZqYGrBIVt19N8ycCaeeCi0p3lfr2GNhueWq9wWfdFADFsmiQvpday34znfSPdaIEXDwwXDddTBrVrrHajJqwCJZ9Oc/w2OP+fCzAQPSP96JJ/rd1H7+8/SP1UTUgEWy6Oc/9xktfvCD2hxvjTX8Tmn//d8wb15tjtkE1IBFsubZZ/3875FH+rnZWjn+ePjsM79KTqpCDVgkayZP9qFhhx5a2+NutBHssYc34IULa3vsBqUGLJIl8+bB1Vf7sLNVV6398U84AebP9y/kpN/UgEWy5Je/9NMAxx0X5/g77gibbw4XXKBp7atADVgkKxYtgssug69+1U8HxGDmzf+553y6I+kXNWCRrLj5Zp9k86ij4tbx7W/76Y+LL45bRwNQAxbJiosv9suOd901bh2DBvmFGXfcAa++GreWjFMDFsmC6dNh6lQfeparg1/bww7zOjQkrV/q4G9SRHp08cU+c/H3vhe7ErfGGj578uWXw6efxq4ms9SARerdu+/6bMUHHOCzXtSLI4+E99+HG26IXUlmqQGL1LsrrvChZ4cfHruSzrbbDsaOhUsu0ZC0PlIDFqlnbW0+9GzHHWGTTWJX05mZfyjMnOnnp6XX1IBF6tldd8Hs2XDEEbErKe873/HTIv/1X7ErySQ1YJF6dsklsPrqMHFi7ErKGzrU75J2003wzjuxq8kcNWCRevXSSz7b8SGH1Oaev311+OGweLGPiJBeUQMWqVeXXeZTDR10UOxKlm2jjWDCBPjVr3wae6mYGrBIPfr0U7jySh9rO2pU7Gp6dsQR8MYbcOedsSvJFDVgkXp0/fU+xrZev3zr6mtfgzXX9HPWUjE1YJF6E4I3srFjfaxtFrS0+A3i77sPnn8+djWZoQYsUm8eeQSeeMKvNDOLXU3lfvhDGDhQQ9J6QQ1YpN5cfDEMHw777x+7kt753Odgn33gqqvgww9jV5MJasAi9WTuXLjlFjjwQFh++djV9N5RR8HHH/vsydIjNWCRevLLX0J7e/3d96FS48fDllv6aYj29tjV1D01YJF68Y9/+FjaPfeEddeNXU3fHXWUX0Ryzz2xK6l7asAi9eK663zG4WOPjV1J/3zrWz52+cILY1dS99SARepBCN6wNt0UdtopdjX9M3Cgj+C47z546qnY1dQ1NWCRenDfffDMMz7jcJaGnnXnkENgueXgootiV1LX1IBF6sGFF8LIkbDvvrErqY5VVvEZPK69FubNi11N3VIDFontqafg7rv9suNBg2JXUz3HHQeLFmn6+mVQAxaJ7ec/9zG/WbnvQ6U23NDvY3zJJT42WJaiBiwS0+zZPqnlQQfByivHrqb6Tj4ZFiyAX/86diV1SQ1YJKYLLvAv3Y4/PnYl6dhqK/jnf/b/zkWLYldTd9SARWKZN8+T4f77+60cG9XJJ8OcOXDNNbErqTtqwCKx/OIXPt38qafGriRdu+0GW2wBZ5/tUxdJBzVgkRjmzYNLL4X99vMvqxqZGZx5Jrz6qlJwF2rAIjGcf77f++FHP4pdSW3ssYdScBlqwCK19vbbPjRr330bP/0WFFLwrFm6VWUJNWCRWjvrLB8RcOaZsSuprT328FERZ57pk46KGrBITb34oo98OOQQWG+92NXUlplfdPLmm7pHREINWKSWTjvNb1JzxhmxK4lj++19BuXzzoP33otdTXRqwCK18tBDPt3QiSf6/GnN6txz/dLkH/84diXRqQGL1MKSJX6P3LXXhn/919jVxLXJJj7l0mWX+ezPTUwNWKQWLrnE73p24YUwZEjsauL76U/9lpVHHNHUc8epAYuk7c03/ZzvV74CX/967Grqw/Dh8LOfwaOP+jT2TUoNWCRNIfiIh0WLfKbgRpjtolomTfIv5Y47DubOjV1NFGrAImm69lq480445xxYf/3Y1dSXXA5+8xv/cDr4YP+wajJqwCJpmTMHjj4att3Wn2Vp663noyLuusubcZNRAxZJw+LFsM8+Pvrhyishn49dUf066iifCfqoo3xi0iaiBiyShtNPh0ce8avedOph2XI5uO46WGEF+Na3mmr6IjVgkWq76Sa/1++hh3oKlp6NGgXXXw/PPw8/+EHTDE1TAxappkcfhe9+F7bZxsf8SuV23tnvFXHTTf4viCbQErsAkYbxwgs+C/Caa8Ltt8PgwbEryp4TToBXXvF7Ray1Fhx2WOyKUqUGLFINzz/vCc7Mv9EfMSJ2RdlkBhdf7CNIDj/cv7w8+ODYVaVGpyBE+uvpp/1b/PZ2ePBBfenWXy0tcPPNfv/gQw7xC1galBqwSH/cfbef7zWDBx6AjTeOXVFjGDTI7xw3caIPTzvuOGhri11V1akBi/RFW5tf3bbnnn4xwV//Cv/0T7GraiyFJnzssTB5Muy+u0/n1EDUgEV6a9YsP+Vw+uk+bvWhh2D06NhVNaZ83keTTJni7/PYsXDbbbGrqho1YJFKffKJ39Vs4439PrZXXw033ABDh8aurPEddBDMmOEfdP/v//n54RdfjF1Vv6kBi/Tkgw98WNQ66/h9bL/xDXj2WR/vq7ub1c7GG8O0aXD++fDww758wAH+d5FRasAi5SxZAv/7v/C97/lVWqeeClts4ZcXX3+9TjnEMnCgjxV+8UU45hg/R7zJJrDDDj7d/YIFsSvsFQtNeAs4Wdr48ePD9OnTY5cRT1ubX0jxyCM+muHee+Hvf/fTC/vt58Ohxo2LXaV0NX++30XtiivgpZf8nPEOO8CECf48bhwsv3xqhzezGSGE8X3+eTVggSZowIsWeTqaPx/eeccH+s+eDS+/7I336afh009929VWg1128XONu+6a6i+wVEkI8Le/wa23wh//6NM/gZ8i2mADH6Gy/vrQ2upXKo4a5ROjrryy//328VSSGrBUxfhVVw3TJ05M/0Dd/f9WeL10fQhLP9rb/bmtrfhYssQfixfDZ5/5Y+FCf3z8MXz0kf+5nNGjYcMNYcwYT0tbbum/sDq3m23vvef35Zg5078wfeEF/7BdtGjpbVtaYNgwfwwZAsst50PgBg3yUx4DBvg2+XzxsfrqcMEFasBSHeMHDgzTazVVenfNrfB66Xqzzo9czp9Lfxnyef8lGTCg+Iuz3HL+GDbMTyMMH+6PVVeFkSO98Y4e7dtKc2hv93/9vP66P8+b5/8qWrDAP6Q/+sj/FbRwYfGDfNEi/2BfsqT4gd/e7rNb33tvvxuw7gUhbtNNoZFPQYjkcn7qYdSo2JV00CgIEZFI1IBFRCLROWABwMw+Al6IXUcPRgDvxi6iB1moEbJRZxZq3DCEMKyvP6xzwFLwQn++TKgFM5uuGqsjC3Vmpcb+/LxOQYiIRKIGLCISiRqwFEyJXUAFVGP1ZKHOhq9RX8KJiESiBCwiEokasIhIJGrATc7MdjOzF8zsZTM7JXY9BWa2ppk9YGbPmtkzZnZM8vrKZvYnM3speV6pDmrNm9lMM7szWV7HzKYl7+mNZjYwcn3DzexmM3vezJ4zs63r7X00s+OSv+enzewGMxtcD++jmf3GzOaZ2dMlr5V978z9Z1Lvk2bW4/1L1YCbmJnlgUuArwIbA/uaWb1M67sEOCGEsDGwFXBEUtspwP0hhPWB+5Pl2I4BnitZ/hlwYQhhPWABcGCUqoouAu4JIWwEbIbXWjfvo5mtARwNjA8hjAHywD7Ux/v438BuXV7r7r37KrB+8jgYuKzHvYcQ9GjSB7A1cG/J8qnAqbHr6qbW24Ev41frjUpeG4VfQBKzrtHJL+HOwJ2A4VdvtZR7jyPUtyLwKskX7iWv1837CKwBvAGsjF8cdiewa728j0Ar8HRP7x3wK2Dfctt191ACbm6F//EL5iSv1RUzawW+AEwDRoYQ3kpWvQ2MjFRWwWTgJKA9WV4FeD+EsCRZjv2ergPMB65MTpNcbmbLU0fvYwhhLnA+8DrwFvABMIP6eh9Ldffe9fr3SQ1Y6pqZDQVuAY4NIXxYui54zIg2jtLM9gTmhRBmxKqhAi3AOOCyEMIXgE/ocrqhDt7HlYCJ+IfF6sDyLP3P/rrU3/dODbi5zQXWLFkenbxWF8xsAN58rwsh/D55+R0zG5WsHwXMi1UfsC3wL2b2GvBb/DTERcBwMyvcZyX2ezoHmBNCmJYs34w35Hp6H3cBXg0hzA8hLAZ+j7+39fQ+luruvev175MacHP7G7B+8m3zQPyLjz9Ergnwb5SBK4DnQggXlKz6AzAp+fMk/NxwFCGEU0MIo0MIrfh7978hhP2BB4BvJpvFrvFt4A0z2zB5aQLwLHX0PuKnHrYysyHJ33uhxrp5H7vo7r37A3BAMhpiK+CDklMV5cU68a5HfTyA3YEXgVeA02PXU1LXdvg/7Z4Enkgeu+PnWO8HXgLuA1aOXWtS747Ancmf1wX+CrwM3AQMilzb5sD05L28DVip3t5H4CzgeeBp4BpgUD28j8AN+Hnpxfi/Jg7s7r3Dv4C9JPldegof1bHM/ffrUmQz2w3/J1ceuDyEcF6fdyYi0mT63ICTMaQv4kOD5uD/nN03hPBs9coTEWlc/bkh+5eAl0MIswDM7Lf4N5ndNuARI0aE1tbWfhxS+uv9933y1zXX7Pz6jBkz3g0hrFpY/nLuW73/ZO4627Hluix2s77L61bYTy7Xeb/JcnF9YRblLvvJ5TuWO7bN5zvvq+P1ws/6c8h1OXa+cw2h4/XulpPn5OdC3sqvT57bW7r+XOGZzsvJYdq7rO+6XNiuuL7zftpbOq9f6rlwnK7btYRu9hs6r0+eSV6nJWDJn63FR+rl88lzsjxggI80G5BvA2Bgiz8Pyheeff1yLYsBGJw8L5/3KeaXyyfLLZ8BMCTnrw/L/wOAocnzCrmFyfrPkmV/fVjHs+9nmIVk2d+EobnBAORWe6mb6bz7rj8NuNyYty27bmRmB+NXhbDWWmsxXTPvRnXSSXDxxUtPgGxms+NUJNK8Up+SKIQwheSemePHj9e9LyNbvBgGDEhp54XTWYV0GZJrE5KEGtqTJJTrsr69c4ItnBaz9mR9IWUmy4XUaYVLH3Jd9kNb8pzvSHTWlrxWSMIFbe2dFi0ZGBTo/HohCRdqKlzDZHRdLugoLllPl/XJ2uQyg/alfhMLW3b+yVyy3N7NcleFd6Q92S6XbNdedusydXVTT3G/5Y8byv7Zf6qNrqrchvq7uyQJ0164BiRJ0P3cbdlD9eNn63oMqZS3aBEMjHprGBEp6M9nRccYUrzx7gPsV5WqJDWffgrLL5/yQeoqCeeTbZMSlISTZSXhHnVJwmkk4D6XGEJYYmZHAvfi/5f/JoTwTNUqk1R8+CEMHRq7ChGBfn5GhBDuAu6qUi1SA++9B6usUqOD1UUSLp4P9m2TEpSEk+X0k3DXL34ym4RToEuRm8zcubDaarGrEBGowSgIqR9tbfDaa/DNb/a4aXXFTMJlRkb4tkkJSsLJspJwDErATeTpp2HJEhgzJnYlIgJ1+Zkgafnzn/15hx0iFRAjCS9jjLBvm5SgJJwsKwnXkhJwE7n2Wth446UvQxaROOros0DSNHUq/O1vfhlydDVMwpVcLefbJiUoCSfL1UvClYwR7rzcPElYCbgJLFoEhx0GI0fCAQfErkZECurgM0DSFAKcfjo88QTcdhusUOnlPGbFpJqWGiTh3tw3ApSE00jCvblarvNy4ydhJeAGFoLf/ez88z0BT5wYuyIRKaUE3KA++ACOPhquvhqOOAL+8z/7sJOOZJrhJNyHO6j59kkJSsLJcn+ScO/vG9F5uXGTsBJwA7rtNh/tcO21cMYZ/sVbTn/TInVHCbiBPPYYnHsu3HEHbLqpN+IvfrGPO7NcSRLNcBLux72EffukBCXhZLkvSbjvd1DrvNx4SVi5KOPa2uDWW2HbbWGbbeAvf4FzzvEZL/rcfEWkJpSAM+q55+DGG+G66+Dll6G11c/zfv/7VbzdZGGutQwn4WrMquHbJyUoCSfLlSfhatxLuPNy4yRhNeAMefFFb7q/+53f18HMLys++2z4xjegRX+bIpmiX9k69vHH8PDDcN998Kc/wVNPedPdbjv/Ym2vvWDUqHSObTnrSJpZTsLVnF/Ot09KUBJOlntOwtWcVaPzcvaTsBpwHVm8GKZN84Z7//1++fCSJT6H27bbwuTJfivJNdaIXamIVIMacETz5nnDLTweeww++cQD3xZbwAknwC67ePNdbrna11dIlplOwinMtOzbJyUoCSfLTZCEU1D/FTaIzz7zy4GnTvVmO3UqvPqqr8vnYexYmDQJJkyAHXeElVeOWq6I1IAacAo+/hiefNIb7syZ/vzkk35THPBTCFtt5ZcHb7mlp93UZyrurZJxwJlOwqnMtAxKwqVHVxLuq/qtLCPeeadzo505E156qdgrVl4ZvvAFOOYYb7ZbbgmjR8etWUTqgxpwhRYuhGef9ZEIpY+33y5u09rqzXb//WHzzf3Po0cXA1um5KyY+rKchFOZaRmUhLvsv6MaJeHeqL+KImtrg1mzlm60L7/c8bvJ4MF+r4Vddy022s02g+HD49YuItnS1A143rylG+0zz8Cnn/p6M/j85/0Lsn328eexY2G99ZYOOA2pcC41w0k4lZmWS/ejJNx5/x3VFJNwf+4l3HmfjZeE66eSFH36qTfWrs123rziNquu6s31oIOKjXaTTerwyzERaRgN14D//nf/Iuzxx4vPL75YDEuDB/u07HvsUWy0Y8f6dD1SZFa8iizLSTiVmZZBSbgXSbgWMy2XVrz0cv0m4fgV9FEI8OabSzfb118vbrPmmjBuXOfTB5//fJOcPhCRupeZBtze7qcRHnrI74/w0EPw1lvF9RtsAFtv7bM/fOEL/hgxIl69mZfLFVNXlpNwCjMt+/rkmErCnZa7ypHOTMvLOmaWknDdNuDFi2HGjGKzfeQRWLDA162xhl8tttVWnnA32wyGDYtarohIr9VVAw7Bbyh+9dV+y8UPP/TXN9jAb7e4/fZ++8XW1oyOrc0Ss470l+kknMJMy6Ak3Jsk3Jcxwp1q6qaWRkjCPR7RzNYErgZG4v8NU0IIF5nZysCNQCvwGrB3CGFBX4qYNQuuucYb76xZPvJgr73ga1/zWy+utlpf9ioiUt8qaflLgBNCCI+b2TBghpn9CfgecH8I4TwzOwU4BTi5twX8x3/AiSd6WJkwAc4809Ouhn9FVjKjcKaTcBozLYOSsJJwVfR4pBDCW8BbyZ8/MrPngDWAicCOyWZXAQ/SywY8ZYo33732ggsv9FELIiLNolet3sxagS8A04CRSXMGeBs/RVGxV1+FQw+F8ePh+uv9puNSP8ysI91lOgmnMdMyKAkrCVdFrudNnJkNBW4Bjg0hfFi6Lvj/nWV/i8zsYDObbmbT58+f3/H66NGw004+fveuu/pWvIhIllXU4s1sAN58rwsh/D55+R0zGxVCeMvMRgHzyv1sCGEKMAVg/PjxHU16wAC47Tb48pf9FMSuu/oNySdO9KvVJLKcdaQ5JWGUhPuRhNOYabnzdtlNwj0mYDMz4ArguRDCBSWr/gBMSv48Cbi9twcfNgzuuQdOPtnvzbDPPj7i4eCDfdxve09/EyIiGVZJa98W+C7wlJk9kbx2GnAe8DszOxCYDezdlwKGD4dzzoGf/hQefBCuugquuw5+/WtYaSUfhrbDDj4GeNw4T85SA5brSG9KwigJ9yMJp3Ev4U41dVNL9ZNw9VUyCuIvLP33UTChWoXk8z4MbcIEuPRSuP12eOABvxLujjt8myFD/HLjwgUZX/qShquJSHbV1ZVwBUOH+qwS++/vy2+/7VfIFe4DcdZZxSGeG27o930YN654DwhNaFkFpTNiKAkrCRfW9iEJpzmrRqeauqmlnpNwXTbgrlZbDb75TX8AvP8+PPoo/O1vPoriL3+BG24obr/22sWGPG6cz1qx+uq6fFlE6ksmGnBXw4fD7rv7o+Ddd70ZFx6PP+6jLApBZ6WVirek3HRTfx4zRjfx6VYuT8dnv5KwknDnvfYqCddifrlONXVTS3+TcBoy2YDLGTHCh7R9+cvF1z76CP7v//zx1FM+NfzVV/vrBa2tnW/MPnas3/yn6b/syxmFZqBGrEbcv0bc+4s1ikeup0ZcfQ3TgMsZNsxHUWy3XfG1EGD27KWnJ7r7bliS/MUNHAgbbbR0Y87sDMciUpcaugGXY+apt7XV77ZW8Nln8MILnpILTfnPf/YhcQXDh/tpi66nMlZYodb/FenzS5ELS0rCSsL9ScJ9v2y5eOT4STgNTdeAuzNokDfUTTft/PqCBfD0053T8nXXFe9VDD5LcmF6+sLzaqspLYvIsqkB92CllXzc8fbbF18LAd54w9PyE0/44/HH4eabi9t87nPFhrz55rDFFt6oM9OU8/mOtKUkjJJwv5Jw/2/gUzxyYyVhNeA+MIO11vLHnnsWX//gA2/KM2d6U545Ey64wKdXAh+fvOWW/thqK7+QZKWV4vw3iEh8asBVtOKKS6flRYvg2Wdh+nSYNg2mTvX7XxQC1gYbeDPeckvYZhs/BZLLld9/TZl1pCslYZSEO36geZNwGtSAUzZwYPE0xA9/6K99+KE35KlTvSnfc48PjwMfTrfTTrDLLn5Z9rrrZui0hYj0ihpwBCusADvv7A8oDo176CG4/35/3HSTr2ttLd4jY7fdanjKorTrKwkrCffPc1UpAAAO+klEQVQjCac55b2vz24SVgOuA6VD4w44wH/XX3ih2IxvuQWuuMIvDvnKV2Dvvf2+ySuuGLtyEekPNeA6ZOYXgmy0ERxxBLS1+X0vbrkFfvc7+OMf/dTGbrt5M95rrxRuYp/PLZWelIRREu5DEk5zeqPikbOZhOvh6x7pQT7vX9T94hfw2mvw2GPemGfMgO98x28+9NOfwnvvxa5URHpDDThjzLwZX3ABvP463HefjzE+4wyfVfrII/18cr/lcp6G8jk/aOkjn/dxwmaYmae6nPkNfHL54rLl/JEsF7fP+aOwv2S5Y31HDV32U3gPclZMmP5Cp/Ud+01TCJ2Tdmgvpl08CXck9dL17cEfHbsJnobb2/1R2G+yXFyfPLrup70tefhyx/Ztbf4o7K/waGv3R7J/aw9Ye8nxC+uT7a293dNwW4C2csvJo609eQQsWddpfXsgt8QfxZ8pPEgeyXK7J/9cWyBXsr7rcmG74np/FPaTW+Ipt7j/Lo/Ccbput8T80WW/aVADzrBczr+cu+uu4pROU6bAxhvD5Mn++yci9UsNuEGMGQO/+Q289BLsuCMcd5yPK37mmb7tL+RKkmqWk3DalIR7TMKlKTjLSTgNasANZu214c474frrYdYsb8KPPBK7KhEpR6MgGpAZ7Luv34Zzl1186Nof/uCnKyqWy3V8Q25dP6ezNDqiFiMjSvev0RGJ4uiIvt1Bzfdauhx/dET1KQE3sDXX9Is71l3Xp3OaOzd2RSJSSgm4wY0cCbfe6veYOOggH0Nc0anRfPGzOdNJuJZjhEv3ryScaKdv940o3bJeknD1KQE3gfXWg3PO8Vk/HnggdjUiUqAG3CQOPdRvh3nppRX+gFlx9EMyOiHkLHujI2KMEQaNjigdAVHBGOEsjI5Igxpwkxg8GCZNgttvh08+iV2NiIAacFP5yld84tGpU3veNpQm0Awn4ZIFJeFISbg3V8vVcxJOgxpwE9l6a3/+61/j1iEiTqMgmsiKK/qoiFdeqWDjvBW/0S58E5/v/HmdhdERdXEHtdL9N+HoiHRmWl76J2s1OqKalICbTGur38RHROJTAm4yI0bAW2/1vF3I5ZYe25nFJFxP9xIu3X8TJeF05pcr3TK7SVgJuMmstBIsWBC7ChGBXiRgM8sD04G5IYQ9zWwd4LfAKsAM4LshhEXplCnVMmQILFxYwYal54CznITrcVaN0v03QRJOd6bl0i2zl4R7k4CPAZ4rWf4ZcGEIYT1gAXBgNQuTdAweXGEDFpHUVZSAzWw0sAdwNnC8+SDMnYH9kk2uAs4ELkuhRqmiAQNg8eKetws5K8kVGU7C9Ty/XOn+GzgJpzHTMtQ+Caeh0gQ8GTiJ4ju0CvB+CKFwgd4cYI1yP2hmB5vZdDObPn/+/H4VK/3X0uIXY4hIfD0mYDPbE5gXQphhZjv29gAhhCnAFIDx48enHCekJ7lcMTAuS8jnoCPJJq9lMAlnYqbl0v03YhJOYaZlqH0STkMlpyC2Bf7FzHYHBgMrABcBw82sJUnBowHdbTYDKm3AIpK+HhtwCOFU4FSAJAGfGELY38xuAr6Jj4SYBNyeYp1SJWaVhbmQN0ozAGQzCVdlVg1QElYSTkV/9n0y/oXcy/g54SuqU5KkqdIGLCLp69WVcCGEB4EHkz/PAr5U/ZIkTRXfACxvJWe+MpyEqzm/HCgJ9yEJV3N+OV9Pl/XJ2pSTcBp0JZyISCS6F0STqTQBdx4HXJDBJJzGTMugJNyLJJzGTMu+ni7rk7UZSsJKwCIikSgBS1mhZDaJLCfhVGZaLlmvJJyUoSTcJ0rAIiKRKAFLWe0tttRMsFlMwtWYVQOUhPuVhKswq0bn5YLaJuE0KAGLiESiBCxlhZx1JIZMJ+Eqzi8HSsLNnITToAQsIhKJErCU5feCcJlOwinMtAxKwr1KwinMtNx5uSDdJJwGJWARkUiUgKWskIeu2SGbSTiFmZZBSbgXSTiNmZZ9+9om4TQoAYuIRKIELGX5OeDyCSBbSTiFmZZL96MknKzvPgmnMdOyr691Eq4+JWARkUiUgKUsTy3LPheWhSScykzLoCTcqyRc/ZmWS2uqXRKuPiVgEZFIlIClrPa8lYx/zG4STmOmZVAS7l0STmF+OYiQhKtPCVhEJBIlYCkr5MvdDSqDSTiFmZZLa1ISriAJpznTMmQ6CSsBi4hEogQsZZWeA1YSVhLuVxJOYaZlqH0SToMSsIhIJErAUla5c8CZTMJpzLQMSsK9SMJpzC8HtU/CaVACFhGJRAlYygolH81KwigJ9yMJpznTMmQ7CSsBi4hEogQsZYX80q9lMQmnMdOyHzM5hpKwPy8rCacxvxzUPAmnQQlYRCQSJWApqz3f/adzlpJwKjMtg5Jwb5JwmjMtQ82ScBqUgEVEIqkoAZvZcOByYAweUH4AvADcCLQCrwF7hxAWpFKl1FzIG+1JFlUSBiXhvifhVGZahghJuPoqTcAXAfeEEDYCNgOeA04B7g8hrA/cnyyLiEiFekzAZrYisAPwPYAQwiJgkZlNBHZMNrsKeBA4OY0ipfbaWyCXZIBMJ+FUZloGJeHKk3AaMy1DjCRcfZUk4HWA+cCVZjbTzC43s+WBkSGEt5Jt3gZGlvthMzvYzKab2fT58+dXp2oRkQZQyTngFmAccFQIYZqZXUSX0w0hhGBmZT+KQwhTgCkA48ePT/njWqrF7wXhspyE05lpGZSEe5GEU5hp2bdPSqhREk5DJXueA8wJIUxLlm/GG/I7ZjYKIHmel06JIiKNqccEHEJ428zeMLMNQwgvABOAZ5PHJOC85Pn2VCuVmiq9Ei7LSTiNmZZBSbhXSTiVmZah1kk4DZVeiHEUcJ2ZDQRmAd/H/+//nZkdCMwG9k6nRBGRxlRRAw4hPAGML7NqQnXLkXpR/l4QLktJOI2ZlkFJuFdJOJWZlqHWSTgNuhJORCQS3QtCygrL+GjOUhJOZX45UBIuqCAJpzLTcul+apSE06AELCISiRKwlNXe0vNssFlIwlWdVQOUhLuzjCScykzLUPMknAYlYBGRSJSApaxOV8JlOgn34r4RoCTcX+WScAozLfv65Ji1SsIpUAIWEYlECVjKCi0BOhKty2YSrvy+EaAkXDWlSTiFmZahMZKwErCISCRKwFKWXwnXOXlmMQmnMdOyb68kXJEQenffCKjbJJwGJWARkUiUgKWszrMiZzcJV2VWDVAS7o80ZlqGmifhNCgBi4hEogQsZYV8KEmgHa8mz0rCdGyvJFwxJeGlKAGLiESiBCxl+Thgl+UknMpMy6Ak3B9ZTcIpUAIWEYlECVjKKk3ABVlMwmnMtAxKwlWhJKwELCISixKwlJcPdJd9lISVhKsqK0k4BUrAIiKRKAFLeSXngLOdhPtzL+FidUrCSsJpUAOWsqzMKYhsNuLqTW/Ueb0acWrqtBGnQacgREQiUQKWsqylncLnc7aTcN9vZVn6U0rCSsJpUAIWEYlECVjKyufbSz73s5uE05zyvvN6JeHU1E0Srj4lYBGRSJSApax8SzGJZTkJpzvlfbE6JeEmSMIpUAIWEYmkogRsZscBP8QDwFPA94FRwG+BVYAZwHdDCItSqlNqbMCAJXT93yOLSTjdKe9Lj64k7MsNnIRT0GMCNrM1gKOB8SGEMfjf9D7Az4ALQwjrAQuAA9MrU0Sk8VR6DrgFWM7MFgNDgLeAnYH9kvVXAWcCl1W7QIljQL507GN2k3CaU96X/pSScOMn4TT0uOcQwlzgfOB1vPF+gJ9yeD+EUPjfcg6wRrmfN7ODzWy6mU2fP39+daoWEWkAPSZgM1sJmAisA7wP3ATsVukBQghTgCkA48ePT/HjUKppYEu5q3+UhJWE+5CE00zBXmTxWJ2OXd0knIZKsvUuwKshhPkhhMXA74FtgeFmVviNHA3MTalGEZGGVMk54NeBrcxsCLAQmABMBx4AvomPhJgE3J5WkVJ7g/LLuv5dSVhJuBdJuBbng0v3n1YSTkEl54CnATcDj+ND0HL4KYWTgePN7GV8KNoVqVUpItKAKhoFEUL4MfDjLi/PAr5U9YqkLgzK9xA7gSwk4d7cN6J0WUm4ikm4liMjSvdf7SScAl0JJyISie4FIWUt17K4F1vXbxJOc6JPUBKuKAnHGCNcuv8qJeE0KAGLiESiBCxlDe5VAi6ovyRciynvQUl4WUk46tVypfvvZxJOgxKwiEgkSsBS1vL5/tzYrn6ScF/uoObrlYT9uMlx+pGE6+K+EaX772sSToESsIhIJErAUtZy+b6cA+5KSVhJuDgOOPNJOAVKwCIikSgBS1nLt3xWxb0pCTd1Eq7HewmX7r/SJJwCJWARkUiUgKWsIblFKfzfoSTclEm4nmfVKN1/D0k4DUrAIiKRKAFLWcPy/yguZDgJpzHTsq9XEvbjJsdZRhLOxPxypfvvJgmnQQlYRCQSJWApa2hpAi7IYBLuz72EO++zu2MqCftxk+OUScKZmmm5dP9dk3AKlIBFRCJRApayVsgt7H5lppJw/2fV6LzP7o6pJOzHTY5TmoTTmGnZd0SquibhFCgBi4hEogQsZQ3JVXAlXCaScPXmlyvuc1nHVBL24ybHyeXSmWkZap+EU6AELCISiRKwlLVCrswoiO7UcRLuz9VynbZTEvbt+pCEU5lpuWR9zZJwCpSARUQiUQKWsob1JgEX1GESTmOmZV9WEobKknAqMy37C53WZzEJKwGLiESiBCxlDcv1Y0YMJeEyx23iJJzCTMu+mP0krAQsIhKJErCUNcwC9CcFg5Jw2eM2YRJOYaZlaIwkrAQsIhKJErCUNSzXAu1J3MpwEk5zVo1O2ykJ+3Zlk3AKMy1DQyRhJWARkUiUgKWsobnBQDIWWElYSbhfSTiFmZZL95PhJKwE3GS22gqOOSZ2FSICYKGGnwZmNh+YXbMDSm+sHUJYNXYRIs2kpg1YRESKdApCRCQSNWARkUjUgEVEIlEDFhGJRA1YRCQSNWARkUjUgEVEIlEDFhGJRA1YRCQSNWARkUjUgEVEIlEDFhGJRA1YRCQSNWARkUjUgEVEIlEDFhGJRA1YRCQSNWARkUjUgEVEIlEDFhGJRA1YRCQSNWARkUj+P4VunBEI1yzqAAAAAElFTkSuQmCC\n",
      "text/plain": [
       "<Figure size 360x360 with 3 Axes>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "#%% plot the distributions\n",
    "\n",
    "pl.figure(1, figsize=(6.4, 3))\n",
    "pl.plot(x, a, 'b', label='Source distribution')\n",
    "pl.plot(x, b, 'r', label='Target distribution')\n",
    "pl.legend()\n",
    "\n",
    "# plot distributions and loss matrix\n",
    "\n",
    "pl.figure(2, figsize=(5, 5))\n",
    "ot.plot.plot1D_mat(a, b, M, 'Cost matrix M')"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Solve Unbalanced Sinkhorn\n",
    "--------------\n",
    "\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 5,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "It.  |Err         \n",
      "-------------------\n",
      "    0|1.838786e+00|\n",
      "   10|1.242379e-01|\n",
      "   20|2.581314e-03|\n",
      "   30|5.674552e-05|\n",
      "   40|1.252959e-06|\n",
      "   50|2.768136e-08|\n",
      "   60|6.116090e-10|\n"
     ]
    },
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAWAAAAFgCAYAAACFYaNMAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADl0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uIDIuMi4yLCBodHRwOi8vbWF0cGxvdGxpYi5vcmcvhp/UCwAAIABJREFUeJzt3XmcXGWV//HP6S2dDUJIDCEBGoZtgKBglE34iUFBkcERRRAlKrLLJiqCzojjiBsCygAaAQdk+WFwAEWWAURRlkhC+LGvgUCAkIBJIGtvz++P89zqrttddKe3p6r6+3696nVzb92qOlXQp0+feu7zWAgBEREZejWpAxARGa6UgEVEElECFhFJRAlYRCQRJWARkUSUgEVEElECFikzZvYLM/u3QXjeW81sZi/P/bOZfXl975P1owQsFcPMgpltnTt2tpld1Wl/nJldYmaLzWy1mT1qZl/sdP/KTrd2M1vTaf+IQY7/C2b2t57OCyEcF0L4Xh9f4ywzeyG+n0Vmdl2n5/1oCOGKvjyvDI661AGIDBQzawDuBJYAewCLgBnAFWa2UQjhvBDCmE7nvwh8OYRwZ4p4u2NmtSGEtj4+dibweWC/EMLzZrYJ8C8DGuAAMrO6EEJr6jhSUgUs1eTzwObAp0MIL4QQWkIItwEnA/9hZhus7xOa2X+b2cXxz/eVZnavmW1iZheY2TIze8rMdul0/jfN7Hkze9vMnjCzf43H/xn4BbBHfJ7lnZ7/EjO7xcxWAfvGY/8Z7z/DzOaYWV3cP97MHjezxm7CfR9wewjheYAQwuIQwqxOsRVaB1k1bmbnxvfxgpl9tMRnMNnMHjGzr3c6vEX8LN42s/81swmdzv+XGOPy+Jr/3Om+F+N7egRYZWZ18djX4musMLPrSry/qqMELNXkw8CtIYRVueO/AxrxqrgvDgW+DUwA1gH3Aw/F/euB8zqd+zywN7Ah8F3gKjObHEJ4EjgOuD+EMCaEMK7TYz4LfB8YC+RbFD+Jr/ltM9sGOAf4XAhhbTdxPgAcaWZfN7PpZlbbw/vaDXg6vo8fA5eZmXU+wcy2BP4C/FcI4Se5mL8IvAtoAL4Wz98WuBY4FZgI3AL8If51kjkcOBAY16kCPhQ4ANgS2Bn4Qg+xVwUlYKkmE4DX8gfjD/kb8f6+uCGEMC8mvRuAtSGEK2Or4DqgUAGHEGaHEF4NIbSHEK4DngXe38Pz3xRCuDc+piixhhDagSPxKv73wI9DCPO7e5IQwlXAScD+eNJcYmZnvMPrLgwh/Cq+jyuAycCkTvfvANwNfKdzJR39OoTwTAhhDfBb4D3x+GeAP4YQ7gghtADnAiOBPTs99uchhJfjYzsfezWE8A/gD52er6opAUslaQPqc8fqgZb47zfwJFIk/vk+Id7fF693+veabvY795WPNLOH45/fy4Gd6Dnxv/xOd4YQXsQTYRNwUQ/nXh1C2A8Yh1fc3zOz/UucvrjT41bHf47pdP8RwCt4lV/yscDqTo/bFFjY6Xnb8fc3pdP53b3fUs9X1ZSApZK8hCehzrak4wf+TuCjZjY6d84h+J/xDwxmcGa2BfAr4CvAxrHN8BiQ/VlfaurBd5yS0MwOxNsnd+EtiR7F/vds4BH8l0BfnI3/0rqmF+2MzKvAFtlObGlshifyQnh9jKfqKAFLJbkO74VONbMaM9sPOIiOCu03+MiH2WbWZGb1sfr7OXB2CGHFIMc3Gk8uSwHi8LfOye91YGquH/qO4pdblwJfBmYCB5nZx0qc+wUzO9DMxsbP56PAjsCcPr0b/8vi0/j7utLMepMvfgscaGYzzKweOB3/5XdfH2OoakrAUkn+A/9B/huwDP/i6IgQwmMAIYR1wH74n7hzgLfwL8i+lfsCaVCEEJ4Afop/Sfc6MA24t9MpfwIeBxabWW/bIbPwHvEtIYQ3gaOAS81s427OfQs4C/9LYTn++RwfQuhx7HEpIYRm4JN4b/jynpJwCOFp4HPAhXj1fBBwUHweyTFNyC4ikoYqYBGRRJSARUQSUQIWEUlECVhEJBFNxiMATJgwITQ1NaUOQ6SizJs3740QwsS+Pl4JWABoampi7ty5qcMQqShmtrDns0pTC0JEJBElYJHhZMUKuO8+WLgQdA1AckrAItUuBJg9G3baCcaNg732gqYmeNe74JvfhLffTh3hsKUELFLNFi+GGTPg0EOhtha+/3246Sa4+GLYd1/40Y9gu+3g9ttTRzos6Us4kWq1aJEn30WL4KKL4NhjPQlnjj8e5syBo4+Ggw6Ca6+FQw5JF+8wpApYpBq98grsvbdXwP/7v3DCCcXJN7PbbvDXv8L73udV8uzZQx/rMKYELFJtWlo8mb7xBtx1l/d838mGG3oLYo89YOZMeOyxoYlTlIBFqs43vuEjHS67DKZP791jxoyB66/3ZHzIIfDWW4MbowBKwCLV5eab4YIL4OSTvQpeH5tsAtddB88/DyeeODjxSRElYJFqsXKl93p32gl+0sf55/fZB846C666Cu64Y2Djky6UgEWqxdlnw8svwy9/CQ29XvWoq7POgm228WS+dm3P50ufKQGLVINHHvHWw9FHw5579nz+O2ls9HHCzz0HP/jBwMQn3VICFqkG3/gGbLAB/PCHA/N8++0Hn/mMtzJefXVgnlO6UAIWqXR33+3DyM46C8aPH7jnPeccaG2F//iPgXtOKaIELFLJQvD5HKZOha98ZWCfe6ut/Oq5Sy+FZ54Z2OcWQAlYpLLdcAP8/e/w3e9673agffvb/rz/9m8D/9yiBCxSsUKA//xPH7Fw5JGD8xqTJsGpp/olyk8+OTivMYwpAYtUqltvhfnz4cwzoW4Q59U69VQYOXLgvuCTAiVgkUqUVb+bbw6f+9zgvtaECXDMMXD11bBgweC+1jCjBCxSif7yF7j/fh9+Vl8/+K/3ta/5bGo//vHgv9YwogQsUol+/GNf0eJLXxqa15syxWdK++//hiVLhuY1hwElYJFK88QT3v/9yle8NztUvvpVWLfOr5KTAaEELFJpLrjAh4Ydd9zQvu7228OBB3oCXrNmaF+7SikBi1SSJUvgyit92NnEiUP/+qefDkuX+hdy0m9KwCKV5Be/8DbAaaelef0PfhDe8x447zwtaz8AlIBFKkVzM1xyCXz0o94OSMHMk/+TT/pyR9IvSsAileL6632RzZNOShvHZz7j7Y8LL0wbRxVQAhapFBde6Jcd779/2jhGjPALM/7wB3jhhbSxVDglYJFKMHcuPPCADz2rKYMf2+OP9zg0JK1fyuC/pIj06MILfeXiL3whdSRuyhRfPfnSS2H16tTRVCwlYJFy98YbvlrxkUf6qhfl4itfgeXL4dprU0dSsZSARcrdZZf50LMTTkgdSbEPfACmTYOLLtKQtD5SAhYpZ21tPvTsgx+EHXdMHU0xM/+lMH++96dlvSkBi5SzW26BhQvhxBNTR9K9z33O2yL/9V+pI6lISsAi5eyii2DTTeHgg1NH0r0xY3yWtNmz4fXXU0dTcZSARcrVs8/6asfHHjs0c/721QknQEuLj4iQ9aIELFKuLrnElxo6+ujUkbyz7beHGTPgl7/0Zeyl15SARcrR6tXw61/7WNvJk1NH07MTT4SXX4abb04dSUVRAhYpR9dc42Nsy/XLt7yDDoLNNvOetfSaErBIuQnBE9m0aT7WthLU1fkE8XfeCU89lTqaiqEELFJu7r0XHn7YrzQzSx1N7335y9DQoCFp60EJWKTcXHghjBsHRxyROpL18653wWGHwRVXwFtvpY6mIigBi5STV16B3/0OjjoKRo9OHc36O+kkWLnSV0+WHikBi5STX/wC2tvLb96H3po+HXbbzdsQ7e2poyl7SsAi5WLtWh9L+/GPw1ZbpY6m7046yS8iue221JGUPSVgkXJx9dW+4vCpp6aOpH8+/Wkfu3z++akjKXtKwCLlIARPWDvvDPvumzqa/mlo8BEcd94Jjz6aOpqypgQsUg7uvBMef9xXHK6koWelHHssjBwJP/tZ6kjKmhKwSDk4/3yYNAkOPzx1JANj4419BY+rroIlS1JHU7aUgEVSe/RRuPVWv+x4xIjU0Qyc006D5mYtX/8OlIBFUvvxj33Mb6XM+9Bb223n8xhfdJGPDZYulIBFUlq40Be1PPpoGD8+dTQD74wzYNky+NWvUkdSlpSARVI67zz/0u2rX00dyeDYfXf4P//H32dzc+poyo4SsEgqS5Z4ZXjEET6VY7U64wxYtAh+85vUkZQdJWCRVH7yE19u/swzU0cyuA44AN77Xvj+933pIilQAhZJYckSuPhi+Oxn/cuqamYGZ58NL7ygKjhHCVgkhXPP9bkfvv3t1JEMjQMPVBXcDSVgkaG2eLEPzTr88OqvfjNZFbxggaaq7EQJWGSoffe7PiLg7LNTRzK0DjzQR0WcfbYvOipKwCJD6plnfOTDscfC1lunjmZomflFJ6++qjkiIiVgkaF01lk+Sc2//3vqSNLYe29fQfmHP4Q330wdTXJKwCJD5Z57fLmhr33N108brn7wA780+TvfSR1JckrAIkOhtdXnyN1iC/j611NHk9aOO/qSS5dc4qs/D2NKwCJD4aKLfNaz88+HUaNSR5Pe977nU1aeeOKwXjtOCVhksL36qvd8P/IR+MQnUkdTHsaNgx/9CO67z5exH6aUgEUGUwg+4qG52VcKrobVLgbKzJn+pdxpp8Err6SOJgklYJHBdNVVcPPNcM45sM02qaMpLzU1cPnl/svpmGP8l9UwowQsMlgWLYKTT4a99vKtdLX11j4q4pZbPBkPM0rAIoOhpQUOO8xHP/z611Bbmzqi8nXSSb4S9Ekn+cKkw4gSsMhg+Na34N57/ao3tR7eWU0NXH01bLABfPrTw2r5IiVgkYE2e7bP9XvccV4FS88mT4ZrroGnnoIvfWnYDE1TAhYZSPfdB5//POy5p4/5ld770Id8rojZs/0viGGgLnUAIlXj6ad9FeDNNoObboLGxtQRVZ7TT4fnn/e5IjbfHI4/PnVEg0oJWGQgPPWUV3Bm/o3+hAmpI6pMZnDhhT6C5IQT/MvLY45JHdWgUQtCpL8ee8y/xW9vhz//WV+69VddHVx/vc8ffOyxfgFLlVICFumPW2/1fq8Z3H037LBD6oiqw4gRPnPcwQf78LTTToO2ttRRDTglYJG+aGvzq9s+/nG/mODvf4d//ufUUVWXLAmfeipccAF87GO+nFMVUQIWWV8LFnjL4Vvf8nGr99wDU6emjqo61db6aJJZs/xznjYNbrwxdVQDRglYpLdWrfJZzXbYweexvfJKuPZaGDMmdWTV7+ijYd48/0X3r//q/eFnnkkdVb8pAYv0ZMUKHxa15ZY+j+0nPwlPPOHjfTW72dDZYQeYMwfOPRf++lffP/JI/29RoZSARbrT2gp/+hN84Qt+ldaZZ8J73+uXF19zjVoOqTQ0+FjhZ56BU07xHvGOO8I++/hy98uWpY5wvVgYhlPASVfTp08Pc+fOTR1GOm1tfiHFvff6aIbbb4d//MPbC5/9rA+H2nXX1FFK3tKlPovaZZfBs896z3iffWDGDN/uuiuMHj1oL29m80II0/v8eCVggWGQgJubvTpauhRef90H+i9cCM8954n3scdg9Wo/d5NNYL/9vNe4//6D+gMsAyQEePBBuOEG+OMfffkn8BbRttv6CJVttoGmJr9ScfJkXxh1/Hj/79vHVpISsAyI6RMnhrkHHzz4L1Tq/7fseOf7Q+h6a2/3bVtbx6211W8tLbBund/WrPHbypXw9tv+7+5MnQrbbQc77eTV0m67+Q+seruV7c03fV6O+fP9C9Onn/Zfts3NXc+tq4OxY/02ahSMHOlD4EaM8JZHfb2fU1vbcdt0UzjvPCVgGRjTGxrC3KFaKr1UcsuOd77frPhWU+Pbzj8MtbX+Q1Jf3/GDM3Kk38aO9TbCuHF+mzgRJk3yxDt1qp8rw0N7u//189JLvl2yxP8qWrbMf0m//bb/FbRmTccv8uZm/8Xe2trxC7+93Ve3vv32fidgzQUhbuedoZpbECI1Nd56mDw5dSQFGgUhIpKIErCISCLqAQsAZvY28HTqOHowAXgjdRA9qIQYoTLirIQYtwshjO3rg9UDlszT/fkyYSiY2VzFODAqIc5KibE/j1cLQkQkESVgEZFElIAlMyt1AL2gGAdOJcRZ9THqSzgRkURUAYuIJKIELCKSiBLwMGdmB5jZ02b2nJl9M3U8GTPbzMzuNrMnzOxxMzslHh9vZneY2bNxu1EZxFprZvPN7Oa4v6WZzYmf6XVm1pA4vnFmdr2ZPWVmT5rZHuX2OZrZafG/82Nmdq2ZNZbD52hml5vZEjN7rNOxbj87cz+P8T5iZj3OX6oEPIyZWS1wEfBRYAfgcDMrl2V9W4HTQwg7ALsDJ8bYvgncFULYBrgr7qd2CvBkp/0fAeeHELYGlgFHJYmqw8+A20II2wPvxmMtm8/RzKYAJwPTQwg7AbXAYZTH5/jfwAG5Y6U+u48C28TbMcAlPT57CEG3YXoD9gBu77R/JnBm6rhKxHoT8GH8ar3J8dhk/AKSlHFNjT+EHwJuBgy/equuu884QXwbAi8Qv3DvdLxsPkdgCvAyMB6/OOxmYP9y+RyBJuCxnj474JfA4d2dV+qmCnh4y/7HzyyKx8qKmTUBuwBzgEkhhNfiXYuBSYnCylwAfANoj/sbA8tDCK1xP/VnuiWwFPh1bJNcamajKaPPMYTwCnAu8BLwGrACmEd5fY6dlfrs1vvnSQlYypqZjQF+B5waQnir833By4xk4yjN7OPAkhDCvFQx9EIdsCtwSQhhF2AVuXZDGXyOGwEH478sNgVG0/XP/rLU389OCXh4ewXYrNP+1HisLJhZPZ58rw4h/E88/LqZTY73TwaWpIoP2Av4FzN7Efi/eBviZ8A4M8vmWUn9mS4CFoUQ5sT96/GEXE6f437ACyGEpSGEFuB/8M+2nD7Hzkp9duv986QEPLw9CGwTv21uwL/4+H3imAD/Rhm4DHgyhHBep7t+D8yM/56J94aTCCGcGUKYGkJowj+7P4UQjgDuBj4VT0sd42LgZTPbLh6aATxBGX2OeOthdzMbFf+7ZzGWzeeYU+qz+z1wZBwNsTuwolOronupGu+6lccN+BjwDPA88K3U8XSK6wP4n3aPAA/H28fwHutdwLPAncD41LHGeD8I3Bz/vRXwd+A5YDYwInFs7wHmxs/yRmCjcvscge8CTwGPAb8BRpTD5whci/elW/C/Jo4q9dnhX8BeFH+WHsVHdbzj8/frUmQzOwD/k6sWuDSE8MM+P5mIyDDT5wQcx5A+gw8NWoT/OXt4COGJgQtPRKR69WdC9vcDz4UQFgCY2f/Fv8ksmYAnTJgQmpqa+vGS0l/Ll/vir5ttVnx83rx5b4QQJmb7H675tGZpEunkjvbZJZbz7rv+JODuxrztlj/JzI7Brwph8803Z65W3k3qG9+ACy/sugCymS1ME5HI8DXooyBCCLNCCNNDCNMnTpzY8wNkULW0QH196ihEBPqXgMt6DKl0r7kZGpJODSMimf4k4LIdQyqlrV4No0enjkJEoB894BBCq5l9BbgdH4Z2eQjh8QGLTAbFW2/BmDGpoxAR6Oey9CGEW4BbBigWGQJvvgkbb5w6Cuk1G/Av3ov14zoA6T9dijzMvPIKbLJJ6ihEBJSAh5W2NnjxRfinf0odyTBk1rdbyeer6dttsOOS9aIEPIw89hi0tsJOO6WORESgnz1gqSx/+Ytv99knbRxVrbfVYXfVaNHdg1Vl1nZ7NLT3shcc4rzzPb1P9ZZ7RRXwMHLVVbDDDl0vQxaRNFQBDxMPPAAPPuiXIcsA6KkCLFHhlqxs8+f3UAFbP/uw+Um4LCuMS1XCWeUbK+iSFXNvKmRVxwWqgIeB5mY4/niYNAmOPDJ1NCKSUQVc5UKAb30LHn4YbrwRNtggdUQVZj0r3S4Vbg+VbaGSrcmdl3/d3P39rYC7VKHt7d3fXzheEw/78S4Vc28q5FLV8TCuiFUBV7EQfPazc8/1Cvjgg1NHJCKdqQKuUitWwMknw5VXwoknws9/njqiCjFQFW883qXCze133G/dPn/H8R7uX0+Wr1y7VMS5+wv7uco4t1+okK3T88XnCO252HvqFw+DylgVcBW68UYf7XDVVfDv/+5fvOX/whWR9FQBV5H774cf/AD+8AfYeWdPxO97X+qoKkSpKqxUxVuiEi1UtLWxSZqvdAvHs8q3eN+y35T5irk2t5/f9lYI3W/zFW5bW9F+aPNq1bKqtS2rgHMVcu5xtLcTQvbZ9LJfPIx6xaqLKlxbG9xwA+y1F+y5J/ztb3DOOb7ihZKvSHlTBVyhnnwSrrsOrr4annsOmpq8z/vFL2q6yfWSr7LWt8ebVbT53m5WsRYq3toSx+N+XawGa3MVcG4/1HZfAYfC63f7LiErKktUwNaW6+1mFW7ct9a2bo9nFW+hQm5v63Lccv3ikFXJsfK1muIRGD32iquoElYCriDPPONJ97e/9XkdzPyy4u9/Hz75SajTf02RiqIf2TK2ciX89a9w551wxx3w6KOedD/wAf9i7ZBDYPLk1FFWmL72envq8ZaqeGNlW/jtmFW6heNZ5Zvtx0o3v63N9mPFW5PbZm+r1KiIrNrMise4zUZDWGt79/tZZdvSVnSc1mw/brOqNu6H1taO423F1XK+ku11RVyFvWEl4DLS0gJz5njCvesuv3y4tdXXcNtrL7jgAvjUp2DKlNSRishAUAJOaMkST7jZ7f77YdUq/wX/3vfC6afDfvt58h05MnW0Vaq3lW+u11vYry2ucEtWvPVxP25DfW3Rtr0hbuv9ddtjpdveUNz7ba+P2/gy2fGQFYlZcZhtc8Whtee3fkJNa3y7rVkFHI+3+Im1LVkl3F503JpjpRsrZFrifmun/VgNd1TFMfh8RZz1jbM3kauIq7E3rAQ8RNat88uBH3jAk+0DD8ALL/h9tbUwbRrMnAkzZsAHPwjjxycNV0SGgBLwIFi5Eh55xBPu/Pm+feQRnxQHvIWw++5+efBuu3m1q5WKB1lvRzv0VPlmlWyXyjdf4eb2G3xbqHRHZPv+em0jctuGrOL1p2mrL97PKuQQXyZkFXG+Es7eX6HnGw/EbdfKN769lrhtDkXbmsJ+rIzXxfcf92vWtsbz4hO1tGJZVRy3XSri1qyKj/30QmUcK99st1AJZ2+mRG+4giphJeB+ev314kQ7fz48+2zH/wPjx8Muu8App3iy3W03mDo1bcwiUh6UgHtpzRp44gkfidD5tnhxxzlNTZ5sjzgC3vMe//fUqVpCK6kePvx+93zrY0ma9Xob6ou3jbHSjRVv2wg/r62xpmjb2uiv0zrC4nnE+7NecNyP20IlXB9HEhR6wll1GPezt5cVhe3x/cSqsibbtvjxmvhXWratbfbjtWtj5bsuvt21Hnfdupri++P7q4mVcc3a1o5quNnLamsuvjowWLYfz8uOx0q40M4ucUVdyVESFVAJKwHntLXBggVdE+1zz3WMPW9s9LkW9t+/I9G++90wblza2EWksgzrBLxkSddE+/jjsHq132/mKwhPmwaHHebbadNg6607CiGpMKXWYsuP8y1V+RbG9+ZGOZSofNsK21j5jvRty6hc5TvKn651ZKyAG7N9r+LaCxVxrPoaYnXXEHuv9W0xvOKRA/liMKsW22Ml3NoSP4+sKo3bmljhZhVv7Ro/v25tbj9u69eEov3arEKuq4G12Vjp3NWC2Tb2hrMxzcSecH4gR9YbzvrctGfPU7mV8LBIwKtXe2LNJ9slSzrOmTjRk+vRR3ck2h131JdjIjJ4qi4B/+Mf/kXYQw91bJ95puOXYGOjL8t+4IEdiXbaNF+uR6rMO/R/S89qlpuLITsvXolm+SvYsv1stMOI4lEOWeXbOjpuR/rztIzKKl/ftsRf9Fkl3BYr37bRsbob6dVf/Ujvo45o9O3oEd6sHVnv+yPrfNsQm7s1VnxVWWtsFq9r9fjWtHrFvqrZt2vWeZN53Vrfb14de91r/HF1q2PluypW6quLe9b1cVxzfRy1EWqtcPVebayAa3JzIPcwdLnnSrgt/4DKqYQrNgGHAK++2jXZvvRSxzmbbQa77lrcPvinf1L7QETKQ8Uk4PZ2byPcc4/Pj3DPPfDaax33b7st7LGHr/6wyy5+mzAhXbxSZkr1frO7c9VyRy+4pnibr5Brs0o4XrFWuLIt9lpz43oLPd+Rucp3dCjato3xsq52jPdHR4325uu4kb4d3+hfVGw8YpUfr/f9DWKTdlRs3tZbcXnYEivg1bFkfavVm83LW7z0fnOdB/SPtb6/fI3fv3qVb1tW5sYzF8YrZxUvcRs/D6OjhM0p/BfJrniLu9mMbYXxvrXF9xdW3chW4bB8L7hyxgmXbQJuaYF58zqS7b33wrJlft+UKX612O67e4X77nfD2LFJwxURWW9llYBD8AnFr7zSp1x86y0/vu22Pt3i3nv79ItNTRpbK4Oky2rExf3KLitRZJVeYb7eWPnm527IX8nWULzNRj20jYqjGkZ75Ts6Vr4TxnilO2nk2wBMGbkcgMkNK/z+Ov9h2bhuJQBja/xxDRT3grMKeFUcVvF2u08ysrTVK5jXWnws5eJ1vnz2a40bArCkwSeZXlbnlXFzbUN8v7H3XfjYch3dUNNpDHK8J6tQs3ko8nMTZ3MQF9aZy62+ke8hFlbfKFEJl7EeE7CZbQZcCUzC++CzQgg/M7PxwHVAE/AicGgIYVlfgliwAH7zG0+8Cxb4yINDDoGDDvKpFzfZpC/PKiJS3npTAbcCp4cQHjKzscA8M7sD+AJwVwjhh2b2TeCbwBnrG8BPfwpf+5oXFDNmwNlne7Wr4V9SlrKrtPLz7ubmYCjMUlY4349ns5i11+W3sVqr9+qtLo7rHRVHOWzYsAaAdzV6BZxVvlMb3gRgSp3XPuNj5TsuVoGNsQ+ajTxoDz5KYlXwnvGKdq+cJ8YKelytHx9T689THyvnmjgGoT2+weUx3JY4nrgtznBW0xZHOhTmmIDW1uxY8bnZjGshm1mtLptrOBt7nVv008CsAAAS2UlEQVR9o1Ahx8dlY7Wzz7jL+ImojHvBPSbgEMJrwGvx32+b2ZPAFOBg4IPxtCuAP7OeCXjWLE++hxwC55/voxZERIaL9eoBm1kTsAswB5gUkzPAYrxF0WsvvADHHQfTp8M11/ik4yJlL1+F5e4uzDrWU3HVw/3ZfAc1cVsXK9psVMOIGq9kGy2OB47bsfG8sTX+oz3S/AerNjcKZFSshBssTvqAV9gtcVzz2uDN6pUN3px+qzX2jEf4dk0cN9yazeq2LhvtEeOOoyNq661TlZ/NZZzN5BYr0WwESXZVXjbmOtu2Fv81UejH58f/ljpexnq9KrKZjQF+B5waQnir833Bx4V0+7+UmR1jZnPNbO7SpUsLx6dOhX339fG7t9zSt+BFRCpZrypgM6vHk+/VIYT/iYdfN7PJIYTXzGwysKS7x4YQZgGzAKZPn15I0vX1cOON8OEPewti//19QvKDD/ar1UT6LYRO/b6s/9d9zVEYW5qNQc1Vutk39vnVhLP1zjpWFc6tKFHY5nqj2SxkhV5prA5b4xwRsQ+6tsV/RFe3eiW7Ko7fXRmHTayujeN5g29HBX/ChmwSXbzCHYFXrO1koyH8/rYe+qBZ77c2V5FnW6uJn1f2seYHQVinOYoLI0qy+6xwjj9X9tj8ZMalBhJnT5T9tyhMHBw35T8aoscK2LzDfRnwZAjhvE53/R6YGf89E7hpfV987Fi47TY44wyfm+Gww3zEwzHH+Ljf9vaen0NEpFL1pgLeC/g88KiZPRyPnQX8EPitmR0FLAQO7UsA48bBOefA974Hf/4zXHEFXH01/OpXsNFGPgxtn318DPCuu3rlLNJfocsKvNnwhPw35bkxqIUVfHPnt2ZrpcV5CmJfs6Yl22YrScRKN86f22V8cFwTri3Otray1ivdN+KIgLrc3A5tsbzMerar63x0xLgaH80wNk7s22hrih7XHB+3Onhl/WabDzta2ubjf7Nxwcvi5BRvt8Teb4u/Tkus0LOK3YoXsejYhk798MJfD0jUm1EQf6PkxYTMGKhAamt9GNqMGXDxxXDTTXD33X4l3B/+4OeMGuWXG2cXZLz//RquJiKVq6yuhMuMGeOrShxxhO8vXuxXyGXzQHz3ux3tve2283kfdt21Yw4ILWgpJZXqBbcX9zILveCsx5utTxa32aq/hTkjCnNDxF5vXTafbvE44I5tTfF+rh/aGn80m+OBf8Txts2x8lwVe8LLG/1KtiWNXrFOqPcr4Tas9Yp3bNxmoyUybbH7uDaW3ivavNL9R2ushJv9yrela337xhrfrohzQ6xbE0v2tR5PNv9vTXy/cZAGNS2hUx/cj2X98sK2vXi/y/jcUn3q9tw44QpUlgk4b5NN4FOf8hvA8uVw333w4IM+iuJvf4Nrr+04f4stOhLyrrv6qhWbbqrLl0WkvFREAs4bNw4+9jG/Zd54w5NxdnvoIR9lkf3y3Gijjikpd97ZtzvtpEl8ql72P0CJ375desHZqgqW6wVnvd/4TX7IVm0oXI1VXGHW9Pa3fSyBCy9XuJIszprWnG39vGVr/Ud21Rrvyb7R6BXrokafw2GDEX4F25h6nw1tdK33gOvisItsNENbrKxbYy97VVxsbmXs9a5sjqMr1o4oer2WrPJd5XHUrYorX6zKVsjwu2uz7VqoXVe8knLWD7fW9qJtVslaW9Zvz+aCyFXGPY1qCOU/+iFTkQm4OxMm+JC2D3+449jbb8P/+39+e/RRXxr+yiv9eKapqXhi9mnTfPIffdlXZfKJuFQrIvvhzRaAzCb/jnfn02r+eP7+mlzSKPy5Ha9OqCksA1/8p3vWuqhdGy/fXZOb0H2kP355oyfGFSO8hVDXECdub/Anrotf3tXWZNvskuIYVvyF0xq3LS3xdeLwt7Z18cvGbImiNdkinDG+1cWJty5bmigu61W/JhSO1a7NlrKPEwQ1Z8vVx19usa1T+GVXYptPyKGnxBzKt0VRNQm4O2PH+iiKD3yg41gIsHBh1+WJbr0VWuMPQ0MDbL9918SsFY5FZCBVdQLujplXvU1NPttaZt06ePppr5KzpPyXv/iQuMy4cd62yLcyNthgqN+F9FmJSjgUpjKMx/NfypWqhLPzsso2d0FH9no1hS+esj+74/Hm+CVWc/alXaxI4+KWrdmSPyOLJ3TPJkRvyy2BlC1Tv6Y+m8w8hlH7zn+OZ62PbEHMrBKvz/azijxbpDPbZot0xmF1dfntmvaOyndta3yuuF2bW64+Ls6ZVcKhUBFnLYpC2V68zfS20i2DSXgywy4BlzJihCfUnXcuPr5sGTz2WHG1fPXVHXMVg6+SnC1Pn2032UTVsoi8MyXgHmy0kY873nvvjmMhwMsve7X88MN+e+ghuP76jnPe9a6OhPye98B73+uJWkm5TPRYCWdfBMW7s8tbCw/PTSJeOJ59kRSHqbUVD1sLsbea9T1r4qQ2tbHXWhsr2rrGeEFGbkmjthFxm03knlsSqLA0UPzJzg9/63yJMNDNZOm+7TydJHQaVtacfaFG3Mb9dblt/MKtdk1b4eKUrOK15h4q33ic9lzvt1QvONf7LXz5Vsa934wScB+Yweab++3jH+84vmKFJ+X58z0pz58P553nyyuBj0/ebTe/7b67X0iy0UZp3oOIpKcEPIA23LBrtdzcDE88AXPnwpw58MADPv9F9st72209Ge+2G+y5p7dAamq6f34ZBD1UwpnCxC5txY8rVL7ZEKpskc78ZD3xG15r84o3q/pC7Ida7OHWZBc2rIlL/cTFPdviJcqFxT7zSx1lUz1mUz9mUz5mPeAs0Pz/W6F4W6iAs0uL24ovoqjJTzLUXFzx1hSGmsW/ANa1Fap9yyrerLdbqHzjZ5Edz74N723l29thZ2XU+80oAQ+yhoaONsSXv+zH3nrLE/IDD3hSvu02Hx4HPpxu331hv/38suyttlLbQqRaKQEnsMEG8KEP+Q06hsbdcw/cdZffZs/2+5qaOubIOOAAtSwGTalxwtnd+d5w3A8UX5LcpTdcG6u2bPKe7BLmOPG5ZT3h5mx5+7i0T32sgLNLmuNy99my94XFP+uzxUCLK+EuUz/mllAKuUq4YxKd4glzOirg3OXD2UUU2f1ZlZtts4srWlo7Kt6ski2MciiugDsq3twoh/WtfPO93zKsfDNKwGWg89C4I4/0/1+efrojGf/ud3DZZX5xyEc+Aoce6vMmb7hh6shFpD+UgMuQmV8Isv32cOKJXgA8+KAn4t/+Fv74R29tHHCAJ+NDDtEk9gOm1KXLXXrDsdIL2SQ6xZODF5ZUb8sm6ckm9cn6ndnCk9mlb7EijtNQZl8EZJVyVvEWlmQvHC+eDKjUfsi/nxKX9Fn+yr3CKIncxPNZBVyYjrOt+LzWTtVuVsHmKtrC+N58r7cwEVL2nN2Paqjkyjejr3sqQG2tf1H3k5/Aiy/C/fd7Yp43Dz73OZ986HvfgzffTB2piKwPC0P4W2L69Olh7ty5Q/Z61a693edM/ulP/VLqkSPhS1+Cr3/dk/L6MLN5IYTp2f6Haz5d/uXDUCr1TajlR0vk1tepKZ60pzDEJbefLUBZeFyu4rXscfltbfHzhdx+l21ND9/ollp6qcsSTcUT5uSXBSqqdvP35Xu8ofg5+9zrLRwfnP9172ifPeBfh6sCrmA1Nf7l3C23dCzpNGsW7LADXHBBx196IlKe1AOuEjvtBJdfDt/5DpxwApx2ms+RfPnlsOOOqaOrAvmqqofREl16xIXHdV8ZF3rF+crYcqMa8hVyqUo7X7HnB5eXqujz7zNXpZaqTrtMFZkdbw9dertdKt3Ca/Wx4i0VewVQBVxlttgCbr4ZrrkGFizwizvuvTd1VCLSHVXAVcgMDj/cp+Hcbz8fuvb733u7QgZIbyvirA1k61cZd3nefO84k6+UC8fzz9PHWitfdebnXchVtZ3uKBzvqdLteKnqr3jzVAFXsc0284s7ttrKl3N65ZXUEYlIZ6qAq9ykSXDDDT7HxNFH+xhiXdo8CEpVY32sjAvnkT8vKlUp5w30xCK5SrfkKKrO1WxPlW6J87reX/kVb54q4GFg663hnHN8qNrdd6eORkQyqoCHieOO84s1Lr64Yw4KGQL9rowzxWMKS1bKhRPi49t6qCr7qodqtdsqtwJXrBhsqoCHicZGmDkTbroJVq1KHY2IgBLwsPKRj/jEUw88kDoSIYQebu3veAttbcW39lB8y98/0LceXq/buHv7GQwjSsDDyB57+Pbvf08bh4g49YCHkQ039FERzz+fOhLpUW8rwRK95OSGWSXbV6qAh5mmJnjppdRRiAioAh52JkyA115LHYUMGFWaFU0V8DCz0UawbFnqKEQE1iMBm1mtmc03s5vj/pZmNsfMnjOz68ysYfDClIEyahSsWZM6ChGB9auATwGe7LT/I+D8EMLWwDLgqIEMTAZHY6MSsEi56FUCNrOpwIHApXHfgA8B18dTrgA+MRgBysCqr4eWltRRiAj0vgK+APgGkI112RhYHkKI60mzCJjS3QPN7Bgzm2tmc5cuXdqvYKX/6uo6VgEXkbR6TMBm9nFgSQhhXl9eIIQwK4QwPYQwfeLEiX15ChlANTVdp24VkTR6MwxtL+BfzOxjQCOwAfAzYJyZ1cUqeCqg2WYrgBKwSPnosQIOIZwZQpgaQmgCDgP+FEI4Argb+FQ8bSZw06BFKQPGTENHRcpFf8YBnwF81cyew3vClw1MSDKYlIBFysd6XQkXQvgz8Of47wXA+wc+JBlMWg1DpHzoSjgRkUSUgIcZVcAi5UMJWEQkESVgEZFElIBFRBJRAhYRSUQJWEQkESVgEZFElIBFRBJRAhYRSUQJWEQkESVgEZFElIBFRBJRAhYRSUQJWEQkESVgEZFElIBFRBJRAhYRSUQJWEQkESVgEZFElIBFRBJRAhYRSUQJWEQkESVgEZFElIBFRBJRAhYRSUQJWEQkESVgEZFElIBFRBJRAhYRSUQJWEQkkV4lYDMbZ2bXm9lTZvakme1hZuPN7A4zezZuNxrsYEVEqklvK+CfAbeFELYH3g08CXwTuCuEsA1wV9wXEZFe6jEBm9mGwD7AZQAhhOYQwnLgYOCKeNoVwCcGK0gRkWrUmwp4S2Ap8Gszm29ml5rZaGBSCOG1eM5iYFJ3DzazY8xsrpnNXbp06cBELSJSBXqTgOuAXYFLQgi7AKvItRtCCAEI3T04hDArhDA9hDB94sSJ/Y1XRKRq9CYBLwIWhRDmxP3r8YT8uplNBojbJYMToohIdeoxAYcQFgMvm9l28dAM4Ang98DMeGwmcNOgRCgiUqXqenneScDVZtYALAC+iCfv35rZUcBC4NDBCVFEpDr1KgGHEB4Gpndz14yBDUdEZPjQlXAiIokoAYuIJKIELCKSiBKwiEgiSsAiIokoAYuIJKIELCKSiBKwiEgiSsAiIokoAYuIJKIELCKSiBKwiEgiSsAiIokoAYuIJKIELCKSiBKwiEgiSsAiIokoAYuIJKIELCKSiBKwiEgiSsAiIokoAYuIJKIELCKSiBKwiEgiSsAiIokoAYuIJKIELCKSiBKwiEgiSsAiIokoAYuIJNKrBGxmp5nZ42b2mJlda2aNZralmc0xs+fM7DozaxjsYEVEqkmPCdjMpgAnA9NDCDsBtcBhwI+A80MIWwPLgKMGM1ARkWrT2xZEHTDSzOqAUcBrwIeA6+P9VwCfGPjwRESqV48JOITwCnAu8BKeeFcA84DlIYTWeNoiYEp3jzezY8xsrpnNXbp06cBELSJSBXrTgtgIOBjYEtgUGA0c0NsXCCHMCiFMDyFMnzhxYp8DFRGpNr1pQewHvBBCWBpCaAH+B9gLGBdbEgBTgVcGKUYRkarUmwT8ErC7mY0yMwNmAE8AdwOfiufMBG4anBBFRKpTb3rAc/Av2x4CHo2PmQWcAXzVzJ4DNgYuG8Q4RUSqTl3Pp0AI4TvAd3KHFwDvH/CIRESGCV0JJyKSiBKwiEgiSsAiIokoAYuIJKIELCKSiBKwiEgiSsAiIokoAYuIJKIELCKSiBKwiEgiSsAiIokoAYuIJKIELCKSiBKwiEgiSsAiIokoAYuIJKIELCKSiBKwiEgiSsAiIokoAYuIJKIELCKSiBKwiEgiSsAiIokoAYuIJKIELCKSiBKwiEgiSsAiIokoAYuIJKIELCKSiBKwiEgiSsAiIokoAYuIJKIEPMzsvjucckrqKEQEwEIIQ/diZkuBhUP2grI+tgghTEwdhMhwMqQJWEREOqgFISKSiBKwiEgiSsAiIokoAYuIJKIELCKSiBKwiEgiSsAiIokoAYuIJKIELCKSiBKwiEgiSsAiIokoAYuIJKIELCKSiBKwiEgiSsAiIokoAYuIJKIELCKSiBKwiEgiSsAiIokoAYuIJKIELCKSiBKwiEgi/x+uBHB0ylWbhAAAAABJRU5ErkJggg==\n",
      "text/plain": [
       "<Figure size 360x360 with 3 Axes>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "# Sinkhorn\n",
    "\n",
    "epsilon = 0.1  # entropy parameter\n",
    "alpha = 1.  # Unbalanced KL relaxation parameter\n",
    "Gs = ot.unbalanced.sinkhorn_unbalanced(a, b, M, epsilon, alpha, verbose=True)\n",
    "\n",
    "pl.figure(4, figsize=(5, 5))\n",
    "ot.plot.plot1D_mat(a, b, Gs, 'UOT matrix Sinkhorn')\n",
    "\n",
    "pl.show()"
   ]
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "Python 3",
   "language": "python",
   "name": "python3"
  },
  "language_info": {
   "codemirror_mode": {
    "name": "ipython",
    "version": 3
   },
   "file_extension": ".py",
   "mimetype": "text/x-python",
   "name": "python",
   "nbconvert_exporter": "python",
   "pygments_lexer": "ipython3",
   "version": "3.6.8"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 0
}
